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ABSTRACT 


NASA and private launch providers have a need to understand and control ignition overpressure 
blast waves that are generated by a solid grain rocket during ignition. Research in accurate com¬ 
putational fluid dynamics prediction of the launch environment is underway. A clearer picture 
is emerging from empirical data which more precisely categorizes all the dissipative mecha¬ 
nisms present in droplet-shock interactions. In this dissertation, water droplets and their effects 
due to vaporization are represented as a control action and two new optimal control problems 
are formulated concerning unsteady shock wave attenuation. A single-phase control problem is 
formulated by representing the effect of droplet vaporization as an energy sink on the right hand 
side of the unsteady Euler Equations in one dimension. Results for the optimal distribution of 
equivalent mass of water vaporized for a given level of attenuation are presented. A two-phase 
control problem consists of solving for the initial optimal water droplet distribution. Results are 
presented for constrained and unconstrained water volume fraction distributions over increasing 
levels of attenuation. New adjoint-based algorithms were constructed which leave the final time 
free and satisfy all first order necessary conditions as well as avoid taking a variation at the 
shock front. 
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CHAPTER 1: 
Introduction 


1.1 Ignition Overpressure 

Ignition overpressure (lOP) is a phenomenon present at the start of an ignition sequence in 
launch vehicles using solid-grain propellants. When the grain is ignited the pressure inside the 
combustion chamber quickly rises several orders of magnitude. This drives hot combustion 
products toward the nozzle and out to the open atmosphere at supersonic speeds. An lOP wave 
is an unsteady shock wave which originates from the exit plane of the nozzle and propagates 
spherically outward near sonic conditions. From previous launch data [l]-[3] it is known that 
overpressures that the body of the rocket experiences are of the order 2:1, or about 15 psi. The 
region below the nozzle, such as the launch pad trench, will experience further compression due 
to displacement of gas along the blast wave’s direction of propagation and overpressures can be 
as high as 10:1, or about 130 psi. The portions of the lOP wave that become incident on the 
rocket body or launch platform components must have an overpressure below a known thresh¬ 
old to avoid costly damage and enable a prolonged lifetime. This is the purpose of the water 
suppression system which is integrated into the launch platform. The current technique used 
by NASA and other launch providers is to spray water into the region around the nozzle before 
ignition. This forces the lOP wave to propagate through water before becoming incident on 
the rocket body or platform components. Through several dissipative mechanisms this causes a 
sufficient decrease to the pressure jump across the shock, preventing damage. 

1.2 Experimental Results on Droplet-Shock Interactions 

A survey of empirical data in the literature was carried out [4]-[21] to understand how droplets 
can be used as a control action in two-phase flow. With blast mitigation as the process of interest, 
an NRL report [7] isolated several desirable and undesirable effects. Drag on the droplets by 
the surrounding gas dissipates momentum. As droplets breakup, the gas does work against 
the surface tension which dissipates energy from the gas. Droplets absorb sensible, latent and 
radiative heat which further transfers energy from the gas. Since water vapor has a higher heat 
capacity than air, after water changes phase from liquid to gas the resulting gas mixture can 
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further absorb heat from the gas. An undesirable effect occurs when liquid turns to water vapor 
so rapidly that the gas mixture density increase dominates the effect on pressure according 
to the constitutive relation. The relative importance of each interaction mechanism isn’t fully 
understood and will change depending on the flow regimes. 

The flow regime of interest will have shocks present with pressure ratios (2:1-10:1) and tem¬ 
peratures up to several thousand Kelvin. Practical droplet sizes have a diameter Dp less than 
500/im. A shock incident on a cloud of water droplets will pass through them so quickly that 
the droplets will not appear to react until after the shock front is spatially isolated downstream. 
Droplets will start to get dragged along in the direction of the gas flow, breakup and vaporize. 
If droplets are large Dp > lOO/rm, then in the presence of strong shocks they catastrophically 
breakup into many small droplets Dp < 25yum. Droplet breakup ceases below a Weber number 
of 12, when inertial forces dominate over surface tension [5]-[8]. Shock attenuation will be 
greatest when the droplets can sink the most energy out of the gas in a given interval of space. 
This will cause the greatest decrease in pressure to the driving gas behind the shock. Pressure 
information will travel at the local sound speed toward the shock front and cause a decrease in 
pressure equal to the diminished driver gas pressure. 

There are general trends in the data which help categorize the relative importance of all of the 
dissipative mechanisms present. First and foremost, droplet vaporization can extract orders 
of magnitude more energy from the gas as droplet breakup can. Secondly, the latent heat of 
vaporization is the most significant dissipative mechanism [8]. Lastly, it has been shown [12] 
that a non-dimensional droplet-flow parameter can predict overpressure attenuation level, shown 
in Figure 1.1. The horizontal axis is a non-dimensional flow parameter. The vertical axis 
is the decrease in overpressure for a shock interacting with droplets divided by the decrease 
in overpressure when no droplets are present, as measured at a fixed location downstream. 
Shock tube experiments were conducted which varied the droplet size and shock strength. The 
data suggests that maximizing the exposed surface area of the droplets will yield the greatest 
enhancement of lOP attenuation. 

1.3 Recent Research on lOP Attenuation 

For many years, the water injection strategy for handling the shuttle’s lOP transient blast has 
been based on order-of-magnitude estimates relating the amount of energy dissipation required 
to the total amount of water used [22]. More than enough water was used and the results were 
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Figure 1.1: Experimental data on lOP attenuation, Jourdan et al. [12] 

suffieient. For the heavy lift vehicles of the future, lOP blast waves will be more substantial 
and require a better understanding of how the water affects the lOP strength. More recently, 
CFD research is underway at NASA [23]-[26] and in the private sector to predict, with greater 
precision, the launch environment during ignition for various solid grain rockets [27, 28]. 

The first work [29] on optimizing one of these precise CFD simulations was a parametric study 
that looked at water arrangement in the nozzle region and how it affected the maximum lOP 
strength. At NASA Huntsville, Cannabal showed in his dissertation [29] and a later publication 
[30] that attenuation is very insensitive to droplet velocity and demonstrated the existence of an 
optimal water injection arrangement. Water cooling the plume near the nozzle has the greatest 
desirable effect of attenuating the transmitted lOP strength. However, an excessive amount of 
water near the nozzle causes obstruction to the blast wave and intensifies pressure. The results 
suggest an optimal arrangement of water exists but there are still an infinite number of possible 
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water distributions even in one dimension. Trial and error, or cost gradient methods based on 
a few discrete inputs are the only option and can yield only coarse notions about continuous 
optimal water distributions. 

1.4 New Contributions of Current Work 

The objective of this work is to develop a computational tool that can directly calculate a droplet 
optimal control for attenuating a range of blast waves to a desired minimal overpressure. As 
shown in Figure 1.2, the water droplets’ control action will be the result of the initial distribu¬ 
tion of the water volume fraction variable ai{x, 0) and the initial droplet diameter. The optimal 
control ai{x, 0) will be the distribution of water droplets which yields the greatest decrease to 
the jump in overpressure of the transmitted shock while not using any more water than neces¬ 
sary. Since the two-phase system and resulting control problem are quite complicated, first a 
single-phase control problem is formulated where the control action is a distributed energy sink 
behind the shock, a simplification for the way droplet vaporization, the dominant dissipative 
mechanism, affects the gas flow. A more sophisticated two-phase model [31] is presented and 
an analogous control problem is formulated in which the control takes the form of the free initial 
distribution of water and evolves dynamically according to the two-phase model. These control 
problems are new formulations and the algorithmic solutions developed in order to satisfy all 
necessary conditions have aspects new to computational results for unsteady shock wave atten¬ 
uation. Given the mathematical framework presented in Chapter 2 and the numerical methods 
given in Chapter 3 of this dissertation, engineering results such as those presented in Chapter 4 
can be obtained for a given incident blast wave boundary condition. 
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Figure 1.2: Physical diagram of simulated interaction. 
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CHAPTER 2: 

Mathematical Theory 

2.1 Euler System in One Dimension 

The mathematical theory of hyperbolic conservation laws has been well established in the 20th 
century [32]-[36]. In order for a PDE to be an accurate physical system the right conserved 
quantities must be identified and an entropy condition must exist. The conservation equations 
dictate that the rate of change of a conserved quantity within a fixed volume equal the flux of 
that quantity across the volume’s surface. The entropy condition selects which conditions are 
physical when discontinuous solutions are present. 

2.1.1 Conservation Equation 

The time-dependent Euler Equations are a system of non-linear hyperbolic conservation laws 
that describe the dynamics of compressible fluids where the effects of viscosity, heat conduction 
and body forces are negligible. In one spatial dimension there are three conserved quantities: 
mass density p, linear momentum pu and total energy pE per unit volume. 

| c /+| f ( C /)=0 ( 2 . 1 ) 

U = {p{x,t),pu{x,t),pE{x,t))'^ (2.2) 


The state vector of conserved quantities U is a function of a time variable which exists in the 
interval [0, T] and a single space variable x which exists in a one dimensional domain E{U) 
is the non-linear flux of conserved quantities. 


F(U) 


( 


pu 

pu^ + P 


\ 


\ u{pE + P) I 


(2.3) 
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For an ideal gas, the eonstitutive relation that relates pressure to internal energy is the ideal gas 
law. 


P = pe (7 - 1) 


(2.4) 


The ratio of specific heat is defined 7 = ^ where Cp is the specific heat of the gas at constant 
pressure, divided by that at constant volume C^. It will be convenient to have the relation be¬ 
tween internal energy and temperature for a calorically perfect gas, e = C^T. A final definition 
relating internal and total energy closes the system from which all thermodynamic quantities 
can be calculated. 


pE 



(2.5) 


2.1.2 Quasilinear Form 

Equation 2.1 can be written in quasi-linear form by creation of the Jacobian matrix A = 


dUi 

dt 


+ Aij{U) 


dU^ 

dx 


0 


(2.6) 


The elements of the Jacobian matrix are given in Equation 2.7. 


MU) 


0 

{pu)^ 7-3 
2 

-7 {pu) {pE) 




+ (7 - 1) 


{pu) 


1 0 ^ 

(7-1) 

i{pE) 3^ i{pu) 


(2.7) 
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2.1.3 Primitive Form 


The Euler system can be written in the primitive variable basis U = (p, u, P)^ by left mul¬ 
tiplying Equation 2.6 by M~^ = and inserting MM~^ between the Jacobian matrix A 
and the term m is the linear gas velocity and P is the gas pressure. Defining the matrix 
A = M~^AM gives the primitive form of the Euler system in Equation 2.8. Details of the 
transformation between bases for a single dimensional fluid flow are given in Merkle [39]. 


dili 

dt 


+ A 






0 


The elements of A are given in Equation 2.9. 

0 ^ 

1/P 

u J 


A{u) = 


u p 
0 u 
0 'jP 


( 2 . 8 ) 


(2.9) 


2.1.4 Characteristic Form 

It can be shown that A has eigenvalues {u,u + c,u — c) where c = y/'jRT is the speed of 
sound in the gas. Any system is hyperbolic if the eigenvalues of the Jacobian matrix are real 
and distinct. Clearly this is true in the case of the Euler System in ID. The specific gas constant 
of air is denoted as R and the gas temperature as T. When the corresponding left eigenvectors 
multiply dU the characteristics equations are obtained. 


dP — pcdu = 0 along 
dP — c^dp = 0 along 
dP + pcdu = 0 along 


dx/dt = u — c 
dx/dt = u 
dx/dt = u + c 


( 2 . 10 ) 


Solving the Riemann Problem for the unsteady Euler Equations amounts to constructing the 
rest of the flow field between these characteristic waves, either approximately or iteratively up 
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to arbitrary precision [40]. 


The method of characteristics takes advantage of the property that hyperbolic systems have 
distinct finite wave speeds. If an observer is traveling at a characteristic speed the observer will 
measure an invariant quantity of the flow. The ID Euler system in characteristic form is given 
in Equation 2.11. It is obtained in a way analogous to the transformation between Equations 
2.6 and 2.8. The matrix A is diagonal and it’s elements (eigenvalues) A* are the real and distinct 
wave speeds {u,u + c,u — c). 


dU ^dU 
dt dx 


0 


( 2 . 11 ) 


Traveling at the characteristic speeds, an observer would measure the Riemann invariants U = 
{S/R,Tc + u,Tc — uY as constant. S is the entropy, R is the specific gas constant, c is the 
speed of sound and T = 2/(7 — 1). Traveling at the speed of a shock wave means that the shock 
will never traverse the observer and hence the observer will never see an increase in entropy. 
Indeed, the Riemann invariants are useful for analyzing isentropic portions of compressible 
flows such as rarefaction waves. 


2.1.5 Well-Posedness 

Eor hyperbolic conservation laws, a well-posed problem requires real and distinct eigenvalues 
of the system’s Jacobian matrix [37]. This is equivalent to the existence of a symmetrizer matrix 
S that has real, positive definite elements and has the property that A A is symmetric. It can be 
shown that the ID Euler System has a symmetrizer matrix as long as density, or temperature 
remain positive [43]. This assumption is built in to the notions of either quantities therefore the 
system is always symmetrizeable and hence the Riemann Problem for the Euler Equations is 
well-posed. 

2.1.6 Entropy and Shocks 

Erom the definition of entropy and the primitive form of the mass and momentum conservation 
equations it can be shown that entropy is constant along dx/dt = u. That is, regions of smooth 
flow are isentropic. The presence of a discontinuity in velocity necessitates that a shock wave 
has formed where the jump in flow quantities are determined uniquely by the Rankine-Hugoniot 
Jump Conditions. Eor a shock wave moving left to right, denote the state just behind the shock 
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Ul and the state immediately upstream of the shoek Ur. The jump quantities, denoted with 
square braekets [ ], are given in Equation 2.12. 


dF {U) ^ [F (U)] = F (Ur) - F {Ul) 
d{U)^[U] = UR-UL 


Rearranging Equation 2.12 gives the shoek speed in terms of the jump across the shock. Equa¬ 
tion 2.13 is the Rankine-Hugoniot jump conditions obtained without any specification to a fluid. 
This condition is a property of any shock wave in a hyperbolic system. 


dx 

Trr shock — 


[F ju)] 
[U] 


(2.13) 


Applying Equation 2.13 to the Euler Equations moving at the shock speed gives the jump in 
primitive variables in a coupled system of algebraic equations which can be solved iteratively 
as will be described in Chapter 3. In addition, an alternative numerical method of Godunov will 
be described which is conservative and, with a stability criteria, obeys Equation 2.13 without an 
iterative solver. 

Since the eigenvalues are increasing functions of the state variables, the the flux is convex with 
respect to U. Eax [32],[33] showed that this requires characteristics on either side of a shock to 
run into that shock. Eor a system with a right moving shock the entropy condition is given in 
Equation 2.14 for k number of waves where 1 < k < n, and n is the number of characteristic 
waves. 


dx 

~g^\shock ^ 

dx 

\{Ur) < -^\shock < \+i{Ur) 


(2.14) 


Intuitively, this agrees with the physical notion that entropy only increases. Eor the ID Euler 
system n = 3. Since a shock requires supersonic conditions the shock speed must be faster 
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than the fastest rate at whieh any information ean travel in the ambient upstream Right state, 
so it must be true that k = 3. This means that the entropy eondition also requires exaetly one 
charaeteristie in the Left state behind the shoek to be faster than the shoek itself. The shock 
is driven by the upstream pressurized gas which must communicate its presence via the u + c 
characteristic in the Left state. This is not hard to believe since the jump conditions show that 
static temperature increases across the shock which raises the speed of sound in the Left state. 
Figure 2.1 shows the instantaneous wave pattern of a shock wave formed from a discontinuous 
and stationary initial condition. 



Figure 2.1: Instantaneous wave patterns of a Left and Right state and shock wave which obey 
the entropy condition of the Euler system in ID 

From the first law of thermodynamics the entropy between states ’L’ and ’R’ is given in Equation 
2.15. 


Sl — Sr = —R ■ In 


Pol 

Pqr 


(2.15) 


Pol and Por are the total pressures to the left and right of the shock respectively. Erom manip- 
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ulation of Rankine Hugoniot Jump conditions it can be shown that for supersonic flows there is 
necessarily a decrease in total pressure, verifying by Equation 2.15 that indeed, Lax’s mathe¬ 
matic entropy condition on the Euler system in one spatial dimension yields the thermodynamic 
principle that entropy can only increase across a shock. 

2.1.7 The Riemann Problem 

The Riemann Problem is an initial value problem where the initial data has a discontinuity. 


U{x,0) 


Ul X < 0 
Ur a; > 0 


(2.16) 


Eor gas dynamics, the Riemann Problem is the general form of the shock tube problem [40]. 
Two stationary gases, one at significantly higher pressure and possibly temperature, are sepa¬ 
rated by a thin wall. As the thin wall ruptures, a shock wave propagates into the gas held at 
lower pressure while a rarefaction wave travels through the gas held at higher pressure. Initially 
a density discontinuity exists between the two gases and this contact surface propagates into 
the lower pressure gas at the speed determined by the jump in velocity at the shock. The shock 
travels at the speed of the gas behind it plus the speed of sound of the low pressure gas. The 
rarefaction wave travels in the opposite direction of the gas at a speed equal to the speed of 
the gas minus the speed of sound in the high pressure gas. The wave pattern of the shock tube 
problem at the instant the diaphragm bursts is shown in Eigure 2.2. 

It is worth pointing out that the fastest wave speed may heu — c since the speed of sound may be 
much greater in the high pressure gas. In Chapter 3, the maximum wave speed will determine 
the maximum allowable timestep discretization and, therefore, computing expense. 

By algebraic manipulation of the Rankine-Hugoniot Jump Condition and the isentropic equation 
of state for an ideal gas it can be proven [40] that solving the Riemann Problem for the Euler 
system amounts to finding the root p* of Equation 2.17. 


f{p*, Ul, Ur) = hip*, Ul) + fnip*, Ur) + [U] = 0 (2.17) 

The rest of the primitive variables u*, p* can be uniquely calculated from p* and the correct 
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Star Region 



Figure 2.2: Solution regions of the shoek tube problem and wave pattern at the initial instant 
the diaphragm bursts 


definition of fi and /r. The present analysis will be restrieted to a left moving rarefaetion wave 
and a right moving shoek wave. The isentropic relation between pressure and density holds 
across rarefaction waves. Solution to the flow in the U^l region is derived from the constant 
Ul state. By combining the isentropic relation and the Riemann invariants /l for a left moving 
rarefaction can be written in terms of the constant left state Ul and the free variable p* in 
Equation 2.16. cl is the speed of sound in the constant Ul state. 


( 2 . 18 ) 

fn in the U^r region is determined by the constant Ur state and the Rankine-Hugoniot jump 
conditions. Applying Equation 2.13 to the ID Euler Equations in a frame of reference moving 
along with the shock gives the jump conditions and, with some rearranging, yields /r for aright 
moving shock wave. 
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1/2 


(2.19) 


/h(p.C/h)=(p-Ph) 

The constants Ar and Br are given in Equation 2.20. The flow can be uniquely determined at 
any time in the star region given a physical initial condition. 


Ar = 


2 

(7 + 1 ) Pr 


Br = 


( 7 - 1 ) 
(7 + 1 ) 


Pr 


( 2 . 20 ) 


The exact solution requires an iterative procedure at every sample point in space and at any 
desired moment in time. This is not practical for large domains or for the added complexity of 
another phase. The difficulty in accuracy of any approximate solution to the Riemann Problem 
comes from any discontinuities at shock fronts. Here finite difference methods fail and a finite 
volume method based on the integral form of the conservation equation, rewritten in Equation 
2.21, is better suited. 



U{x, t 2 )dx 


r^2 fh ft2 

/ U{x,ti)dx+ / F{U{xi,t))dt — / F{U{x 2 ,t))dt 

' X\ Jti Jti 


( 2 . 21 ) 


No discretization or approximation has been made in Equation 2.21 and it is true for any general 
domain [xi,X 2 \ x [^ 1 ,^ 2 ]. 

Numeric schemes based on primitive variables fail at shock waves. If the jump in quantities 
across the shock are wrong, the shock speed will be wrong and therefore, over time, the lo¬ 
cation will be wrong as well. Only schemes based on conservative variables have acceptable 
accuracy near shocks. A numeric method is conservative only if the approximate numerical flux 
is constructed in such a way that each quantity is conserved at every time step. In particular, 
Eax [36] showed that conservative methods that are also convergent will converge to the cor¬ 
rect weak solution at shocks. With that in mind, Godunov’s Method with a suitable choice of 
numerical flux will be described in Chapter 3 and used for all results shown in Chapter 4. 
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2.2 Two-Phase Flow Model 


2.2.1 Assumptions 

The vast majority of the two-phase flow model is based on the work of Saurel and Abgrall 
[31],[41]. It eonsists of six balanee laws and a seventh non-eonservative PDE for the volume 
fraetion variable. The model was ehosen beeause it maintains a hyperbolie strueture and is ap- 
plieable to physieal phenomena of interest, water droplet-shoek interaetions. This is aehieved 
by eonsidering both phases or fluids as eompressible, by assuming many droplets per eell and 
therefore homogenized interfaee eonditions and by solving the same set of equations every¬ 
where in spaee. 

2.2.2 Balance Equations 

The eonservative veetor of the flow is again denoted as U. The subseripts in Equation 2.22 
denote gas or liquid. 


U (oig, ^gPgi ^gPg^gi ^gPg^gi ^iPh ^iPl^h ^iPl^l) 


( 2 . 22 ) 


dag 


+ v. 


dag 

dx 


0 


(2.23) 


Ft 


( 


\ 


agPg ^ 


^gPg'^g ^ 


{ ° 1 


^gPg'^g 


agPgU] + agPg 


Pi 


^gPgPg 

d 

+ 

Ug FgPgEg + OlgPg) 


FV, 

dag 

Oil Pi 

dx 

OllPlUl 


0 

dx 

OilPlUi 


aipiuf + aiPi 


-Pi 


oiipiEi ! 


^ ui {aipiEi + aiPi) j 


V -p^v,) 



+ S{U) 


(2.24) 


The gas phase obeys the ideal gas equation of state while the stiffened gas equation of state is 
used for the liquid. 


Pg^g 


Pg 

7-1 


(2.25) 
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Pi + T^ni 
7z - 1 


(2.26) 


Piei = 


7 = 1.4 is the gas constant of air, 7 ^ = 4.4 is the analogous constant for water [31] and 
TT; = 6 ■ 10® Pa is the stiffening constant that makes large changes in liquid pressure produce 
almost no changes in density. Equations 2.27 and 2.28 are closure relations for the internal 
energy pe to the total energy pE, per unit volume, for each phase. 


Pa ^9 Pa^a “I" 2 ^^^9 


(2.27) 


piEi = piei + -piu^i 


(2.28) 


2.2.3 Interface Quantities 

The volume fraction will propagate at a mean inter-facial velocity which is a center-of-mass 
estimate given in Equation 2.29. 


_ ^aPa^a P ^iPi'^i ^2 29 ) 

^aPa T ^iPi 

Ei and Pi are volume averages of total energy and pressure respectively at the interface shown 
in Equation 2.30. 


Ei OigEg T (XlEl 

Pi ^a^a ^iPi 


(2.30) 


The sum of the volume fractions of each phase will always equal 1 for all a; G fl. 


Olg T Oil - 1 


(2.31) 


When either volume fraction tends toward zero, the ID Euler system is recovered. 
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2.2.4 Primitive Form 


As with the single-phase Euler system, the two-phase system ean be written in non-eonservative 
form. Again let U be the primitive veetor. 


U (^CXg, Pg, Ug, Pg, Pl, Ul, Pl^ 


(2.32) 


In this basis the Jaeobian matrix for the system of Equations 2.23 and 2.24 is given in Equation 
2.33. 


A{u) = 


V 



V 


0 

0 

0 

0 

0 

0 

Ps_ 

“9 

{V^- 

Ug) 

Ug 

Pg 

0 

0 

0 

0 


OgPg 


0 

Ug 

l/Ps 

0 

0 

0 

“9 

Hv^- 

-Ug) 

0 

Pg(Ag 

Ug 

0 

0 

0 

PL 

ai 


Ul) 

0 

0 

0 

Ul 

pi 

0 

9 

Pi-Pi 

»gPg 


0 

0 

0 

0 

Ul 

1/pi 

pi4. 

0^1 


- Ul) 

0 

0 

0 

0 

pvl 

Ul 


J 


(2.33) 


The souree term H{U) whieh multiplies ^ in Equation 2.24 has been eoupled into the Jaeo¬ 
bian matrix and the two-phase quasi-linear system now has the form shown in Equation 2.34. 




(2.34) 


With this fortunate manipulation, the seven eigenvalues can be found by diagonalizing A{U). 
The characteristic speeds are {Vi, Ug, Ug — Cg, Ug 4- Cg, ui, ui — ci, ui -f q), where q is the local 
speed of sound in the liquid phase. The velocity of the phase interface is the seventh charac¬ 
teristic wave speed of the system. All are distinct except locally where they are degenerately 
zero. 


2.2.5 Phase Interactions 

The source vector S{U) is broken up into three separate interactions as defined in Equation 
2.35. 
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S{U) = MH{U) + VR{U) + PR{U) 


( ° 1 


( ° 1 


' a(P,-P,) \ 

m 


0 


0 

rhVi 


Pd 


0 

{Lhv + Pi) + Qi 

+ 

PdV, 

+ 

-pPi {Pg - Pi) 

—rh 


0 


0 

—mV) 


-Pd 


0 

y m {Lhv T Pi) Qi j 


V -PdV. J 


V pPi (Pg - Pi) ) 


(2.35) 


The source vector for mass and heat exchange between the phases is denoted as MH{U). The 
source terms for velocity and pressure equilibration are VR{U) and PR{U) respectively. As 
Saurel and Abgrall point out [31], an interface separating two phases must reach the same 
pressure through microscopic interactions. Indeed, without enforcing this condition, notions of 
thermodynamic properties such as temperature cannot be determined and numerical oscillations 
due to pressure differences will grow without significant artificial dissipation. 

The equilibrium condition Pg = Pi is chosen thereby neglecting the effect of surface tension. 
As stated in Chapter 1, the effect of droplet breakup is omitted in favor of vaporization since the 
latter can be shown to be much more significant of an energy sink to the gas. The microscopic 
pressure equilibration causes a volume and internal energy variation of each phase. Details of 
the solution procedure for the two-phase model are saved for Chapter 3. Isolating the ODE 
^ = PR{U) in Equation 2.35 yields conditions better suited to computation. 


dag 

In 


/i {Pg - Pi) 


(2.36) 


dagPgEg 

dt 


{Pg - Pi) 


(2.37) 


dagpiEi 

dt 


-pPi {Pg - Pi) 


Therefore, Equations 2.36-2.38 simplify to Equations 2.39 and 2.40 


(2.38) 
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(2.39) 


dagPgEg _ 

dt ~ ' dt 


= ( 2 - 40 ) 

dt 

In a closed system, this is just the first law of thermodynamies for an isentropie transformation. 
Time integration of both sides yields eonditions given in Equations 2.41 and 2.42 whose solution 
will be diseussed in Chapter 3. 



{agPgEgf^ - {a g P g Eg)^ = j ^ P^a g (2.41) 

JoO 

(aipiEif^ - {aipiEif = - / ' P^dag (2.42) 

JaO 

The superseript denotes the solution to the ODE for pressure relaxation while ° denotes the 
starting value whieh eomes from the solution to the balanee PDE at the eurrent time step. The 
drag foree Ed exerted by the gas onto a spherieal water droplet, whieh is responsible for the 
VR{U) souree term, is given by the empirieal drag law in Equation 2.43 

Fd = CdpgD^ai {ug -uif (2.43) 

where Dp is the diameter of the water droplet. This is a dynamie variable distributed in spaee 
with the assumption of loeally mono-dispersed droplets. Cd is a eonstant drag eoeffieient. 

The exehange of mass and heat between the phases will be due to vaporization of liquid water 
droplets by the surrounding gas. The rate of gaseous mass produetion in the form of water vapor 
from the liquid phase, m, is defined in Equation 2.44. 


— PH2OV {Pg) Cti 


(2.44) 


The rate of ehange of the volume fraetion of water within a fixed volume of spaee is related 
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to the diameter rate of ehange via Equation 2.45. 


Sp dPp 

Ax^ dt 


(2.45) 


Sp denotes the total surfaee area of droplets per a fixed volume defined in Equation 2.46. 


>S'p = 47r(^) Np 


(2.46) 


Np is the number of droplets per the same fixed volume as is uniquely determined by the volume 
fraetion of water ai for mono-dispersed droplets. 


Np Q-i 





(2.47) 


The evolution of droplet size needed in Equation 2.47 is determined by the rate of vaporization. 
In this work, the Empirieal-Beta Vaporization law is implemented with eonstant values eoming 
from [7]. 


dD^ 

-^ = -HTg) (2.48) 

P {Tg) = 7600 (l + 7.4 ■ 10“^ {Tg - 30077)^) /s (2.49) 

As will be demonstrated in Chapter 4, the range of pressure magnitude for the flows of interest 
to the lOP attenuation problem are from 1-10 atm. Aeeording to steam tables [42], the water 
vapor density will vary by an order of magnitude over this pressure range. Consequently, water 
droplet vaporization will be mueh more effeetive to high pressure shoeks. Equation 2.50 is the 
least-squares quadratie fit to the steam table data for water vapor density as a funetion of the 
surrounding gas pressure. 


pH^OviPg/atm) ^ -.0024344 ■ + .5358 ■ Pg - .077246 kg/m^ (2.50) 
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The latent heat of vaporization of water at its boiling point Lhv = 2.26 MJ/kg. 

The last term in the MH{U) veetor, Qi is the rate of heat exehange between phases at the 
interfaee given in Equation 2.51. 


Q, = hSp{Tg-Ti) (2.51) 

h is the heat transfer eoeffieient. For water droplets of diameter Dp, h = where Nu is the 

Up 

Nusselts Number and A = 0.6 W/mK, is the thermal eonduetivity of water. 

2.2.6 Well-Posedness 

The Jaeobian matrix in Equation 2.33 is symmetrizeable if, in addition to the restrietion in the 
ID Euler System, the volume fraetion of eaeh phase remain non-zero [31]. Sinee our model 
dietates that the same equations be solved in eaeh eell, this restriction is not a problem. 

The eigenvalues of A above are (Vi, Ug,Ug — Cg,Ug + Cg,ui,ui — ci,ui + q). From the defini¬ 
tion for interface velocity and sound speeds it can be seen that all these eigenvalues are real 
and distinct if a few conditions hold. If there are regions in our domain where a single phase 
is completely absent, then the interface velocity will be degenerate to the velocity of the phase 
100% present. Therefore, our calculation restricts the volume fraction of either phase to re¬ 
main above a certain small threshold. As with the single-phase Euler system in one dimension, 
well-posedness can be shown by the existence of a real, positive definite symmetrizer. It was 
shown by Texier [43] that a necessary condition on the positive-definiteness of the symmetrizer 
depends on the positiveness of a certain constitutive variable that should always physically be 
positive such as density, temperature etc. The multi-phase system has the same requirement. In 
addition, compressible flow has been shown to change type from hyperbolic to elliptic where 
the eigenvalues of the system are complex. A famous example is the sonic bubble in aerody¬ 
namics [44]. In that case, the flow is irrotational and allows a velocity potential formulation 
which, in polar coordinates, will change type at the sonic line. Since the flow is steady, a steady 
sonic bubble develops inside in which the flow is elliptic and hence hyperbolic flux-based con¬ 
servative numeric methods would fail. Blast waves, on the other hand, are unsteady and not 
irrotational. Any sonic regions of the flow are dynamic and do not form sonic bubbles. 
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2.3 Optimal Control of Systems of Partial Differential Equa¬ 
tions 

Optimal control theory and the Pontryagin Minimum Prineiple are well-established in mathe- 
matieal and applied engineering literature [45]-[52]. 

2.3.1 Derivation of Necessary Conditions of General Optimal System 

This seetion eontains a general derivation for all the neeessary first-order eonditions of opti¬ 
mality when the dynamieal eonstraints are a system of partial differential equations with free 
or fixed initial and final data, free or fixed spatial boundary data and free final time. In Seetion 

2.4.2 the general eonditions are applied to the unsteady Euler Equations in ID to formulate a 
shoek attenuation problem and to a two-phase flow problem in Seetion 2.4.3 governed by the 
model given in Seetion 2.2. 

To start, a hyperbolie system of balanee laws under a distributed eontrol aetion Zi written in a 
general eonservative-flux form using tensors is shown in Equation 2.52. 

where i = 1,..., m where m is the number of state variables and k = 1,..., n where n is the 
number of spatial dimensions. If the derivatives of the flux veetor is taken with respeet to the 
state veetor U, a Jaeobian matrix for eaeh of the spatial dimensions is defined. 

Atii.U) = (2.53) 

The state dynamies ean then be written in matrix Jaeobian form. It is important to mention here 
that writing the PDE system in Jaeobian form does not imply that this is the form used to solve 
the flow. Eor flows with shoek waves, eonservative methods are neeessary and lend themselves 
to the eonservative-flux form of the equations. However, the adjoint system of equations derived 
below are derived from the Jaeobian form and are always linear in the adjoint variables, making 
their numerieal solution easier. 


23 



S^{U) + z, 


(2.54) 


dUi 

dt 


+ 4 (^) 


dU^ 

dxk 


A non-linear partial differential operator N can be defined so that: 


N{U, z) = ^ + 4(^)^ - S^{U) -z, = 0 (2.55) 

After the dynamical system is known, a cost, or objective, functional must be intelligently 
defined. In the most general framework there can be a cost associated with the initial data J, a 
cost involving the final data K and a running cost L that accumulates over the time interval of 
the control problem [0, T]. So far, initial and final data can still be free or fixed, no restriction has 
been made. The cost functionals describe the objective associated with how these parameters 
are distributed. Let J be the total cost functional. 


J = /(f/(-, 0)) + K{U{-, T)) + r L{U, z)dt (2.56) 

JO 

To incorporate the constraint in the minimization, N is multiplied by a generic continuous vector 
V{x, t) and added to J. The problem then becomes minimizing the augmented cost functional 
J. 


J = /((/(•, 0)) + K{U{-, T)) + r L{U, z) + {V, N{U, z)) dt (2.57) 

Jo 

It is standard to refer to Vi{x,t) as the Lagrange multiplier or adjoint vector. It is a Tri¬ 
dimensional vector since there is an adjoint state for every state variable. 

Let {z*, U* ,V* ,T*) be the optimal control, state vector, adjoint vector and final time respec¬ 
tively. Let a perturbed control function 5zhe added to the optimal control solution so that: 


z = z* -j- eSz (2.58) 

where e> 0 is a small constant. The perturbed optimal control introduces a perturbation to the 
optimal state and adjoint state: 
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u = U* + e6U 


(2.59) 


+ e5V (2.60) 

Another way to see 5U which will be useful shortly is: 

SU = jU [z*+ e5z)\,=o (2.61) 

It has been pointed out in the literature [53]-[55] that if the state variables have shocks, a pertur¬ 
bation is not small in the neighborhood of the shock and does not have the vanishing properties 
as e —)■ 0. A slight increase in the amplitude behind the shock perturbs the speed shock and 
therefore also the location of the shock front. This causes small perturbations to induce varia¬ 
tions on the order of the jump across the shock. The presented method of solution avoids this 

issue, as will be demonstrated later. Since only decreasing the amplitude of a shock wave is 

desired it is apparent that any realistic control action will only slow the shock wave down. The 
target state Q{x) and final time penalty will be constructed in such a way that all variations 
of the solution will occur upstream of the shock front only. Matching the simulated pressure 
profile under control action with the target final state near the shock front will occur by allowing 
the final time to be free. Henceforth, it can be assumed that all variations are taken in smooth 
regions of the flow and that the solution procedure will not depend on a shock location variable 
and adjoint state, nor is a more sophisticated variation required. 

The first order necessary optimal condition that must be satisfied by perturbed control, state, 
adjoint state and final time {6z, 6U, 6V, 6t) is: 

^J{z* + e6z)U=o = 0 (2.62) 

etc 
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= fj {U {z* + e5z)) |,=o + iK {U {z* + efc)) |,=o 

+ ^ /g L {U (^z* + eSz ), z* + e6z) |e=o + iV* + e(5V, N {U [z* + eSz ), z* + eSz)) \e=odt 
= fj {U {z* + eSz)) |,=o + iK {U {z* + eSz)) |,=o 

+ Jo" {U [z* + eSz ), z* + eSz) |f;=o + ^ {V* + e(5V, N {U (^z* + €6z ), z* + eSz)) \f^=Qdt 
+ iL + {V\N{U\z*)))\t=T6t = 0 

= [§l,SU{x,0)) + (§,SU(x,T))+f,^ (i,5^(a:,t)) 

+ (f, dz(x, t)) + {V*,N {U*, l^*)) + {v* + eSV, £N {U {z* + edz ), + edz)) l=^dt 

+ {L + {V%N{U%z*)))\t=Tdt 


(The second term in the second line is zero since N = 0). 

= (f,h[/(a;,0)) + (f,<5[/(a;,T)) 

+ iy*^ lc=o) dt 

+ {L + {V%N{U%z*)))\t=Tdt = 0 


(2.63) 

This is a good place to see structure of the equations. In Equation 2.63 it is easier to see the 
Frechet derivative analogously as a regular chain rule derivative. It is worth noting that in this 
form the equation is valid for any differential operator N, linear or non-linear, PDF or ODE. 
Note that the components of the objectives I, K, L are clearly separated and what remains is to 
carry the perturbation and differentiation with respect to e through on the differential operator 
and then integrate by parts to get the adjoint partial differential system with boundary condi- 
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tions. That term in the above equation can now be written in the Jacobian matrix form with an 
assumption that the eigenvalues are real and distinct. 


{y\ fN {U {z* + e5z ), + efc) |,=o) = 


fe +A{U {z* + eSz)) -S{U {z* + e6z)) - {z* + e6z)] |,=o) 

= IF) + iA A A + |.=o - {V\ §6U) - {V%6z) 


(2.64) 


The last line shows the differential terms separated and makes use of the definition of 6U. The 
last two terms in Equation 2.64 are simple to handle but the first two require more work. It is 
now necessary to use tensor notation. 


Term 1 


A/* d dUi{z*+eSz)\ I 

’ & Ft ) 1^=0 


/T/* 8 dUi{z*+€5z) ) I 

^ dt de ) 1*^=0 

{y;, 


dSUi 

dt 


= \v;ml - (¥.«) 


(2.65) 
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Term 2 


(V*, {U {z* + e5z)) AiA:±A1\ 


dx j 


e=0 


(t?. Uol* + Ai, (u (Z- + £fe)) f/^{z-+ eSz) |„„) 

+ i {A„v-,6Uj) - {^v-,w,) - (aA-^zu, 


( 2 . 66 ) 


vzsu,) + I {A„V-,SUi) - {aZ-^,w, 


Term 3 

[k, p. (U (z- + .fc))) |„„ = (v.-, I^SU,') = (2.67) 

The fourth term simply yields an inner product between the adjoint vector and 5z.li will make 
things clearer to see the two boundary terms in isolation from the differential relations \y*5U^ 
^nAj-^[{AV\5U)]. 

After integrating by parts all the terms in the last inner product, Equation 2.63 can be written, 
for hyperbolic PDEs with source terms, as [56]: 


{§,SU{X,0)) + {%,SU{X,T)) + Ivsuf„+J^ {§,m) + (f ,fe) 

- C-W’A + ((« w - A) K,SU,) - {A,‘§,SUt) - Av-,su,) - iVASz,) 

-I [(AiiV-,SU,)\ \%Z%ldt + (L + (V\N(U\Z-))) \,.TSt = 0 


( 2 . 68 ) 


Grouping the terms that multiply 5U inside the integral gives the adjoint system of PDEs while 
grouping the terms multiplying 5z gives the optimal control condition. Notice the first term in 
the second line multiplies 5Uj not 5Ui but since it is summed over in that term, any index suffice 
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such that it is legal to replace the indices with j so that all terms multiply the same component 
of 5U. Grouping the like variations m5U at the initial and final times gives necessary condition 
on the optimal adjoint variables. 




dV* , / dAij dUk _ 9Aij ^ Y* _ _ dSj t/* j_ 


U_L ( 

dt \ dUk dx 


dx 


O dx dUi 


V* + SU- + 

^ dU-i ’ ' 


(2.69) 


-£ ISg + (i + (C/*,Z-))) |,.T* = 0 


^ fiT \ 

— -V%x,0),SU{x,0)\=0 


(2.70) 


(|^ + V^*(a:,T),(5[/(a;,T)) =0 (2.71) 

Lastly, the final time may be free to vary. If so, a necessary condition for the final time must be 
constructed, known as the transversality condition. The variation of the state at the final time 
will have a first order term expanded in space 6U{x, T) = 6U + St. Therefore, there is 

a change to the final time boundary condition. 


(^6U{x,T) + 


dU{x,T) 

m 


5t, 


dK 

W 







(2.72) 


The first term recovers the same result as for the fixed final time. The second term remains 
and is added to another contribution to the variation of the final time from the integrand of the 
augmented cost functional. This term arises as a result of Leibniz’s rule for a derivative of an 
integral with a moving boundary. When the final time was fixed, the derivative with respect 
to e could be brought inside the integral with no other consideration. Since the final time (the 
boundary of the time integral), is free to vary, Leibniz’s rule dictates that the derivative can be 
brought inside the integral with the addition of boundary terms. For a general function of time 
/ over the time dependent interval[a(t), h{t)] this looks like: 
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M-t) df 

-^■dt + f{b{t))b - /(a(t))d 
it) dt 


(2.73) 


d 

dt 


Mt) 

I a{t) 


f{t)dt = 


In this problem, only the final time is free so d = 0 and the analogous term at b{t) is, 


{L + {y,N{U,z)))\t=T5t (2.74) 

Combining both terms that multiply 5t gives the inner-product: 

(«. (L + (V^-. mu, m I..T + f f + (v^-. f)) = 0 (2,75) 

Recall that TV = 0 meaning that term contributes nothing to the sum at any moment in time. 
The term involving K can be written as a total derivative in time. 






= 0 


The Hamiltonian is defined in general in Equation 2.77. 


(2.76) 


H 


L + 



(2.77) 


Using this definition, the Transversality Condition takes it’s familiar form in Equation 2.78. 


+ =0 (2.78) 

Equations (2.69)-(2.71) and (2.78) gives a set of first-order necessary conditions that an op¬ 
timal system {z*, U*,V*,T*) must satisfy if it can be a candidate for a local minimizer of J 
constrained to N{U, z), the general PDE system under a control action. Initial and boundary 
conditions are either specified or penalized. Eormal sufficient conditions for optimality involve 
going to second-order however it will be obvious by the sequence of iterations of a solution 
procedure if J is being locally minimized or not. 
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2.3.2 Formulation of Single-Phase Optimal Control Shock Attenuation 
Problem 

This single-phase control calculation is meant to give insight into a multi-phase control calcu¬ 
lation where the initial water placement determines the degree of shock attenuation. As stated 
in the Introduction, several interaction mechanisms between the water droplets and the lOP 
wave are present and not fully understood. The dominant dissipative mechanism for the shocks 
of interest is the loss of energy of the gas through vaporization of the water droplets [7]. Ex¬ 
perimental data from droplet-shock interaction in this regime shows that the other dissipative 
mechanisms e.g. drag on droplets, sensible heating of droplets etc. are less significant to lOP 
attenuation. In a multi-phase calculation, the control action would take the form of a liquid 
mass source and the effect of vaporization will take energy out of the gas phase. To most sim¬ 
ply replicate the dominant dissipative mechanism with a single-phase calculation, the control 
will act as an energy sink distributed in space and time behind the moving shock front. No 
sinks or sources in the momentum or mass conservation equations is needed for this simplified 
formulation. 

The cost functional in this context must reflect that a decrease in the maximum jump in pressure 
(the overpressure at the shock front) is most desirable. It should also penalize control action but 
to a lesser magnitude. Therefore, let J, the cost functional be defined in Equation 2.79. 

J= - [ [ z(x,t)'^dxdt -\— [ iPix^T) — Qix))\dx (2.79) 

2 Jo JQ 2 JUs 


The final time T, is not fixed, z{x, t) is the control action, P(a;, T) is the Pressure at the final 
time, Q{x) is the desired final pressure and a and b are weighting constants. The larger b 
is compared to a, the more significant the final time penalty and the less the penalty for using 
control action. Prom the form of Equation 2.56, the individual cost functionals can be discerned. 
There is no free initial data and therefore no need to associate a cost to it, making J = 0 for 
the single-phase formulation. The running cost L{U, z) need only penalize control effort since 
the final time target state sufficiently describes the degree to which a controlled solution had a 
desirable outcome. 
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L{U,z)=^ [ z{x,tY^ 
2 JQ 


(2.80) 


dx 


K{U(;T))=\I {P{X,T) - Q{x))ldx 
I JQs 


(2.81) 


The control action will take the form of an internal energy sink distributed in space and time 
behind the shock front. It only appears on the right hand side of the energy balance equation of 
the ID Euler system. 


U = {p{x,t),pu{x,t),pE{x,t)f 


(2.82) 


d 

(p ] 

d 

^ pu ] 


( 0 \ 

m 

pu 

dx 

pu^ + P 

— 

0 


[pE ) 


\u{pE + P) j 




pE = pe + — 


(2.83) 


(2.84) 


P = pe (7 - 1) 


(2.85) 


Initial conditions are fixed as stationary, ambient air. Stating the density, velocity and pressure 
determines the internal and total energy and the conservative vector. 
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p(a;,0) 

= 1 kg/m^ 

u{x, 0) 

= 0 m/s 

P{x^ 0) 

10^ Pa 

pu{x, 0) 

= 0 kg/m? ■ s 

pe(a;,0) 

= 250000 J 

pE{x,bi) 

= 250000 J 


Therefore, the initial variation in the state veetor in Equation 2.70 is zero. In addition, there is 
no need for an initial penalty / in this formulation so the derivative |^ = 0 and V(x, 0) is left 
undetermined. 

Remark on applying the inlet boundary condition: 

The inlet boundary eondition is explieitly given by the lOP simulation data when the flow is 
supersonie. If the flow is subsonic, non-reflection of the u — c characteristic is imposed [12]. 
In addition, the Monitor Point data was chosen where the flow was nearest to ID; however, the 
2D data did have transverse motion. The data will still give a plausible ID blast wave with the 
inlet boundary condition set in this manner and the goal of the calculation, controlling a range 
of blast waves, can be achieved. At the outlet boundary the flow remains stationary because 
the final time will always be such that the shock wave will not have enough time to propagate 
though the entire domain and reach the outlet. 

In order to determine the optimal control z*{x,t) that minimizes the cost functional J it is 
necessary to define the Hamiltonian of the system and derive necessary conditions using the 
Pontryagin Minimum Principle and the calculus of variations. Recalling Equation 2.77, the 
Hamiltonian for this system is given in Equation 2.87. 


H{U,V,z,t) 




f a ^ dp ^^d{pu) d{pE) 


(2.87) 
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T 

The co-state or adjoint, vector is (Vi, V2, V3) . The flow dynamics will be solved in the con¬ 
servative form of Equation 2.83, yet it is useful to write the system in the quasi-linear form of 
Equation 2.88 to derive the adjoint system. The conservative basis Jacobian matrix was given 
in Equation 2.7. Now the system has a control source vector in the energy balance equation. 


dU, 

dt 


+ ^ijiU) 


du^ 

dx 




( 2 . 88 ) 


To derive necessary conditions as in the general case, the optimal state must be defined 
{U*{x, t), V*{x, t),z*{x, t),T*) and the optimal control perturbed such that = z* + eSz where 
e> 0 is a small constant. The variation in the control will cause variational terms in each of 
the other free variables of the system that must necessarily vanish at an optimal solution. To 
incorporate the constraints of the ID Euler system using the Eagrange multiplier method, each 
conservation law is multiplied by an adjoint variable and added to J, resulting in the augmented 
cost functional J. Then from expanding J in a Taylor series the first order necessary condition 
will be: 


^J{z* + €5 z)U=o = 0 (2.89) 

de 

Grouping the terms of like variational multipliers gives the optimal system. Integrating by 
parts until all derivatives are on the adjoint vector yields the following linear system of partial 
differential equations [11]. 


dt dx y dx 


dAij 



f) f) L 


0 


(2.90) 


Non-trivial elements of the matrix 


dAij 

aUk 


are given in Equations 2.91-2.95. 


d ^ dUk 
dUk^^^' dx 


(7-3) 


V? dp ^ 
p dx 


u d {pu) 
p dx 


(2.91) 
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(2.92) 


d dUk 
dUk' dx 


( 7 - 3 ) 


f u dp 
\p dx 


1 <9(pm) \ 
p dx J 


d ^ dUk 
dUk^"^^' dx 


^ ( 27 E - 3 (7 - 1) m') ^ (3 (7 - 1) 



d {pu) 
dx 


u d {pE) 

p dx 

(2.93) 


d dUk 
dUk^"^^' dx 




3(7-1) 


u d {pu) 
p dx 


I i B(pE) 
p dx 


(2.94) 


—A ■ — = + 7 <9 {pu) 

dUk dx p dx p dx 

dVt denotes the boundary of the spatial domain. The right hand side of Equation 2.90 is zero in 
this formulation because the running cost does not explicitly depend on the state vector. 

Since the final state is not fixed, there is a necessary condition on the adjoint vector at the final 
time. 


v;{x,r)^^K(u(x,T‘)) 


Therefore, in this basis, all 3 derivatives are non-zero in 


(2.96) 
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V{{x,T*) = b{P{x,T*)-Q{x))^-^{x,T*) 

= b{P{x,T*)-Q{x))^- ^^^^P\ ^-l) 


V*{x,T*) = b{P{x,T*)-Q{x))^-J^{x,T*) 

= -b{P{x,T*)-Q{x)),-u{x,T*){^-l) 


(2.97) 


V*{x,T*) = b{P{x,T*)-Q{x))_ 


dP 


d{pE) 

= b{P{x,T*)-Q{x))^-{^-l) 


{x,T*) 


For a free final time, the necessary condition for the optimal final time T* is given by the 
Transversality condition. Define / as a functional to be used in the solution procedure. 


H (U* (x, T *), V* (x, T*), z* (x, T*), T*) + —K (U (x, T*)) = 

at 




/ ^z{x,Ty + b{P{x,T*)-Q{x)) 
Jn 2 




dx = 


= f{T* 


(2.98) 


The necessary condition on the optimal control solution comes from maximizing the Hamilto- 
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nian for an unconstrained control. The unconstrained condition is justified because there is no 
restriction on control magnitude in regions where control is allowed. The integral is true over 
any domain and at any moment in time it should be true point-wise for all t G [0, T]. 

piTT p 

= / az*(x,t) + V^(x,t)dx = 0 (2.99) 

oz Jn 

Equations (2.82)-(2.86),(2.90),(2.97),(2.98) and (2.99) give the complete set of first order nec¬ 
essary conditions for the optimal system. 

2.3.3 Formulation of Two-Phase Optimal Control Shock Attenuation Prob¬ 
lem 

In the two-phase system, the control action takes the form of the free initial data. The first 
state variable ai{x, 0), the initial water volume fraction distribution, will solely determine the 
attenuation at the final time. This is a satisfying formulation as the alternative would involve 
water injection over time as the control. However, water cannot be significantly injected into a 
flow over the interaction time scales of a shock wave. The goal of attenuating the jump in gas 
pressure at the final time is unchanged, so the final time penalty need only adopt the two-phase 
notation. Equation 2.100 gives the cost functional for the two-phase system. 

'^= / ^ {ai{x,0)f dx + j ^{Pg{x,T) - Q{x))\dx (2.100) 

Note that there is no running penalty, thus L = 0 in this formulation. It is then clear what the 
initial penalty and final time penalties are, shown in Equations 2.101 and 102 respectively. 

I{U{;0))= J^^{ai{x,0)fdx (2.101) 

K{U{;T))= [ ^{Pg{x,T)-Q{x))ldx (2.102) 

The multiphase system can be defined in the primitive basis recalled from Equation 2.32. 
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( 2 . 103 ) 


tj (^CXg, Pg, Ug, Pg, Pl, Ul, Pi) 

The state dynamics are of the form of Equation 2.34 with the Jacobian matrix defined in Equa¬ 
tion 2.33. The adjoint system has the form of Equation 2.69. Note the inclusion of the term 

present in the two-phase system but not the single-phase system. This term as well as other 
non-trivial closed form elements of the matrices of the adjoint system are given in Appendix A 
and B for brevity. In the primitive basis, the density, velocity and pressure of the gas and liquid 
have fixed initial conditions which are stationary ambient atmosphere at sea level. 

Pg{x,0) = 1 kg/m^ Pi(a;,0) = 10^ kg/m^ 

Ug{x,0) = 0 m/s ui{x,0) = 0 m/s (2.104) 

Pg{x,0) = 10® Pa Pi{x,0) = 10® Pa 

The first variable, the liquid volume fraction, is not fixed at the initial time and acts as the 
control. The initial cost functional / is then associated with a penalty to the amount of water 
used. 

Eooking back to Equation 2.70 and combing Equation 2.101 and 2.105 gives Equation 2.106. 

<5ff,(a;,0) = | (2.105) 

(2.106) 

After dropping the factor of 6ag(x, 0) the first order necessary optimal condition on V/(x, 0) is 
given in Equation 106. 

f a{l — ag{x,0)) — V/{x,0)dx = 0 (2.107) 

J 

Since this is true over a general domain Vt it must be true point-wise. 

At the final time T, our cost functional K will penalize the gas pressure profile in space if the 
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overpressure at the shock front is greater than that of a predefined target state. K is only an 
explicit function of the fourth state variable U^, the gas pressure. The rest of the state variables 
are free to vary and their variation does not vanish at the final time. 

{ 0 i ^ A 

. ( 2 . 108 ) 

yt 0 i = A 

Using the condition in Equation 2.71 with the definitions in Equations 2.102, 2.108 and 2.109 
gives final time necessary condition in Equation 2.110. 

SP,(x,T)\ (2,109) 

Therefore, as before, the final time optimal condition on V7(a;, T) becomes Equation 2.110: 

(^^ + V:{x,T),6U,{x,T)j =0 ( 2 . 110 ) 

where Q{x) is the predefined target state and 6 is a constant. The optimal condition then be¬ 
comes Equation 2.111: 


[ biPg{x,T)-Q{x)) + V:{x,T)dx = 0 (2.111) 

J 

which must also be true point-wise. 

The purpose of the calculation is to optimize blast wave attenuation. To that end, let the forcing 
function or blast wave condition at a stationary point over time be the left boundary condi¬ 
tion uiit) with a right-ward traveling blast wave. The flow conditions are supersonic and may 
be considered fixed unless reflected shock waves originating from within the simulation do¬ 
main must propagate outward in a non-reflection characteristic boundary condition. The spatial 
boundary term from Equation 2.69 for a linear interval x E [xl,xr] = U is simplified in 
Equation 2.112. 
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i = C? £ {AV* ■ su) dx = . 6U]112 

= A (Ur) V^ -SUr-A (Ul) Vl ■ dill 

= a(u (x^, t)) V*{x^, t) ■ 5U (x^, t)- a(u (x^, t)) V*{x^, t) ■ 5U (x^, t) 

( 2 . 112 ) 

The calculation can be performed in such a way that, after starting with no motion in the initial 
conditions, the blast wave can be simulated as entering at the left boundary, moving right-ward, 
and not reaching the right boundary by the final time. In such a simulation there is no change 
to the state vector at the right boundary for all time. 

{ 0 X = Xr 

0 X = xl supersonic inflow (2.113) 

7 ^ 0 X = Xr subsonic inflow 

From Equation 2.112, the boundary integral from Equation 2.111 reduces to: 
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Because we have all the control effort taking place at the initial time, through the initial data, 
there is no need for control over time. So, z*{x,t) = 0 and since it is fixed, the variation 
6zi{x, t) vanishes as well. 

The final time penalty K can completely communicate the objective of overpressure attenuation 
without the need for a penalty over space and time L for matching a target over space and time. 
All of this means that L{U, z) = 0 and consequently, the first two inner products inside the time 
integral in Equation 2.68 are also zero. The complete set of first order necessary conditions 
for the two-phase optimal system is given by Equations 2.103-2.104, 2.107, 2.111, 2.115 given 
the state dynamics of Equation 2.34 and corresponding adjoint PDE from Equation 2.69. The 
closed form matrix elements for the adjoint PDE are stated in Appendices A and B. 
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CHAPTER 3: 
Numerical Methods 


3.1 Computational Fluid Dynamics for Unsteady Compress¬ 
ible Flow 

3.1.1 Exact Solution to the Riemann Problem 

This section describes a numerical method [40], [57] used for calculating the exact solution 
to the Riemann Problem (Equation 2.16) for the Unsteady Euler Equations of gas dynamics 
(Equations 2.1-2.5). The method takes advantage of the knowledge of the wave patterns shown 
in Eigures 2.1 and 2.2. The subscripts l and ^ refer to the left and right initial data and constant 
states to the left Ul = {pl,ul, Pl) and right Ur = {pR, ur, pr) of the star region respectively. 
Originating from a; = 0 at an initial moment in time, a rarefaction wave moves to the left at 
speed M — c, a contact discontinuity moves to the right at a speed u and a shock wave moves to 
the right at a speed u + c. The speed of sound changes across the star region making it locally 
determined. The exact solution is constructed by sampling the flow at any moment in time t > 0 
on a discrete grid and iteratively solving Equation 3.1 using knowledge of the sample points’ 
location in relation to the wave pattern and given the definitions in Equations 3.2-3.4 previously 
stated in Chapter 2. 

/ {p\ Ul, Ur) = h {p, Ul) + fn {p, Ur) +Ur-Ul (3.1) 



(3.2) 


(3.3) 
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2 


^L = 

Ar = 
Bl = 
Br = 


(7+1)Pl 


(7+l)Pfl 

{^+1)PL 
(7-1) 


(7+1) 


PR 


(3.4) 


For each discrete sample point, the process starts with an initial guess for the pressure p*. Next, 
the value of f{p*, Ur, Ur) is determined based on the known wave patterns. The velocity u* is 
updated according to Equation 3.5 or 3.6 depending on the part of the star region the sampling 
point is located. 


ul = UL-hip*,UL) (3.5) 

U*ji = UR + fR{p*, Ur) (3.6) 

Next, update p* according to Equation 3.7 or 3.8. 

Pi (3.7) 

'>« = (7^) < 3 . 8 ) 

The shock speed in a stationary frame is denoted as S. This solution method is impractical since 
many iterations are required to solve for the root of Equation 3.1 for each discrete point in space 
at each instant in time. The exact solution is shown in red in Eigure 3.1. 

3.1.2 Godunov Method with an Approximate Riemann Solver 

The first-order Godunov Method [58] assumes a piece-wise constant solution in each discrete 
cell volume. The Riemann Problem is then solved locally between each cell at each moment 
in time based on the previous piece-wise solution. This approximate scheme does not require 


44 



an iterative solution and therefore is much less computationally demanding. With the choice 
of a conservative method, the errors introduced by the approximate solution behave in a known 
and acceptable way [59]-[64]. Figure 3.1 compares approximate numerical solutions of the 
shock tube problem to the exact solution, four-tenths of a millisecond in time after the di¬ 
aphragm bursts. The finite difference. Modified MacCormick method [65] introduces dispersive 
error while the conservative Godunov method introduces dissipative error. The dispersive error 
causes the wrong jump condition and therefore the wrong entropy solution, shock speed and 
position as shown in Figure 3.1. Similar results exist in the literature [66]. More oscillations 
in the solution will also cause more problems when equilibrating with another phase. On the 
other hand, dissipative errors are second order in spatial derivatives of U and therefore, are a 
viscous-like error. It is well established that in the limit of vanishing viscosity the correct shock 
solution is achieved. From a physical point of view, a truly inviscid flow is being modeled. 


X 1o' 



Figure 3.1: Comparison of numerical solutions to the Riemann problem .4 ms after diaphragm 
burst in shock tube flow. The initial pressure ratio was 10:1 in each case with a discontinuity at 

X = .5m 


Let the spatial domain x G [xl, xf{\ = be divided into M intervals each equally spaced by 
Ax = {xr — Xl) /M. The locations of the discrete grid points are given in Equation 3.9, for 
j=l,...,M-rl. 
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Xj = (j — l)Aa; 


(3.9) 


The solution will be piece-wise constant between half-cells. As a sensible choice for limits of 
integration, the half-cell locations are define in Equation 3.10, for j = 1,...,M. 


Xj_i/2 = (j - l/2)Aa; (3.10) 

For all cell volumes let the limits of integration be xi = Xj -112 and X 2 = Xj+i/ 2 . Let the time 
domain [0,T] be divided into N intervals each with size At = T/N. The moments in time are 
given in Equation 3.11 for n = 1,...,N. 


r = {n-l)At (3.11) 

For all regions of the discrete space-time domain [a;j_i/ 2 , a;j+i/ 2 ] x the conservation 

law in Equation 2.21 becomes Equation 3.13 without approximation. 




f F{U{xj_i/ 2 ,t))dt - f F{U{xj+i/ 2 ,t))dt 

J ti J ti 

(3.12) 


Discrete approximations to the solution at two discrete moments in time are defined in Equation 
3.13 and 3.14. 


1 r^j+1/2 - 

U{x,F)dx (3.13) 

Ax 

1 /*^? + l/2 — 

U{x,F+^)dx (3.14) 

Ax 

The cell average of the conservative vector is denoted as U{x,t). Locally in each cell, U is 
constant when traveling at the estimated wave speeds. This is a valid assumption as long as 
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the time step is sufficiently small such that no wave is propagated beyond the cell from where 
the speed was estimated to avoid wave interaction which causes inaccuracy. This leads to the 
famous CFL condition given in Equation 3.15. 


At < — (3.15) 

This time interval is limited by the maximum wave speed in each cell S. For t G [t”, t” + At] 
the discrete inter-cell solution and flux are given in Equations 3.16 and 3.17. 


U,.,/2 = U{x,.,f2,n (3.16) 

= (3.17) 

With the discretizations of Equation 3.13, 3.14, 3.16 and 3.17 and the CEE condition on At, the 
integral form of the conservation law in Equation 3.12 gives the numeric form of Godunov’s 
Method in Equation 3.18. 


UF = - F,.y,) (3,18) 

A suitable, conservative approximation to the inter-cell flux T];+i /2 will result exclusively in 
dissipative, not dispersive, numerical error near the shock front. As described earlier, these 
attributes are necessary for accurate jump conditions at the shock. In the limit of vanishing 
viscosity they converge to the correct weak solution. The Rusanov approximation used for the 
inter-cell flux [40], [54] is given in Equation 3.19. 

f;+i/2 ^1{F(U,) + F ((7,+i)) - S, (C/,+1 - U,) (3.19) 

The maximum wave speed in the cell is Sj. For a shock this is given by the Rankine- 
Hugoniot jump condition in Equation 2.13 of Chapter 2. 
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3.1.3 Second Order Godunov Method 


A conservative second order extension to the Godunov method [67]-[70] was implemented to 
solve the compressible flow dynamics in ID under a distributed control action. This method as¬ 
sumes that the solution in each cell is piece-wise linear and projects an intermediate solution on 
a non-uniform grid based on the maximum characteristic speeds from the interpolated solutions 
at each cell interface. Figure 3.2 is a diagram of what the uniform and non-uniform grids might 
look like, based on the wave speed estimates at the inter-cell boundary. 



Figure 3.2: Intermediate non-uniform grid, Kurgarov and Tadmor [67] 


The familiar Godunov integration on the uniform grid is then second order accurate because of 
the intermediate finer grid. This numeric method was chosen because it preserves the conserva¬ 
tive property with only dissipative error and can be adapted to the solution of a ID multiphase 
calculation as presented in Saurel and Abgrall [31]. The essential numerical steps are summa- 
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rized here. 


At a current time step t” the solution vector is known for all Xj . The local spatial derivative 
of the solution vector (f4)” can be calculated by: 



minmod 


uu 


Tjn 

^i+1 


Ax 


Ax 


(3.20) 


The minmod function provides an upwind/downwind switch based on the wave direction of 
motion. 


minmod (a, h) 


- {sign (a) + sign (&)) ■ min (|a|, |6|) 


(3.21) 


The strategy is to first create a finer, intermediate grid to calculate the solution 
at a half time step forward in time with a conservative Godunov method. The second stage is 
to calculate the solution at the next time step on the coarser, uniform grid with a conser¬ 
vative Godunov method. The result is that spatial resolution at the shock front more accurately 
captures the shock front discontinuity while maintaining a conservative second order dissipa¬ 
tive error shape. Other methods can introduce dispersion, losing the conservative property and 
consequently will have oscillations near the shock front which give the wrong jump conditions. 

The upper bound on the maximum wave speed at the inter cell boundary as estimated with a 
first order Taylor expansion of u”. 


<+ 1/2 = max (eig (a (<+ 1 / 2 ) ) , ^9 (<^ 1 / 2 ) ) ) (3-22) 


/■/+ — 
^i+1/2 - ^i+i 


Ax 


{U. 


n 


(3.23) 


^-+V2 = ^ (3.24) 

The superscripts and “ refer to the solution to the right and left of the inter-cell discon¬ 
tinuity. The eigenvalues of A are the eigenvalues of |^, the characteristic wave speed es- 
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timated at the inter-cell boundaries. For the single-phase Euler system the eigenvalues are 

{ug,Ug + Cg,Ug “ Cg) aud for thc two-phasc systcm aTC (yi,Ug,Ug — Cg,Ug + Cg,ui,ui — ci,ui + q). 

The finer, non-uniform intermediate grid is defined by using the bound on the maximum wave 
speed at each inter cell boundary. 


''^j+1/2,1 — "''i-l-1/2 “j-l-1/2 


(3.25) 


— -r” -L . Af 

"''i-|-l/2,r ~ "''i+1/2 “T “j-l-1/2 


(3.26) 


The solution on the non-uniform grid is given in Equations 3.27 and 3.28. 
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i+i/2 




^“i+1/2 


F - F (£/”/«) 


n-rl/2 \ 


(3.27) 


+ 2 flj+1/2) (Ux) 

At/Ax 


1 —At/Ax^ 




['-'j+1/2,1 


F f 7'/'"+1/2 


(3.28) 


Quantities on the uniform grid at the inter-cell boundaries are approximated at the current time 
step in Equation 3.29 and 3.30 and at the half-time step in Equations 3.31 and 3.32. 


£/’Vi/ 2 ,, = u;' + Ax (C/J” (i - 


(3.29) 


C/jVvx, = £/jV. + Ax (£/.)”« (I - 


(3.30) 


Trn+1/2 _ TTu ^^pfrrFi- ) 

^j+i/ 2 ,z - ^i+i/ 2 ,z - [’^j+1/2,1)^ 


(3.31) 
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T tTI-\-1/2 
^■+l/2,r 


ut»n.r-Y^(u 


j+l/2,r 


(3.32) 


The solution is progressed to on uniform grid with a conservative second order update as 
shown in Equation 3.33 and using the definition in Equation 3.34. 


Ax V®J-l/2 + ^j+ 1 / 2 )) + Aa:S'+l/2^i+l/2 + 
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(3.33) 
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1 + — (. 
^ Ax V 


'i+1/2 “i+3/2i 


(3.34) 


3.2 Method of Solution for Two-Phase, Gas-Liquid Coupling 

The two-phase balance equations from Equations 2.22-2.24, the point-wise relations given in 
Equations 2.25-2.28, the interface quantities defined in Equations 2.29-2.31 and source terms 
defined in Equations 2.35 of Chapter 2 are the basis for the overarching numeric method used 
to couple the two-phase system based on the model from Saurel and Abgrall [31]. A survey of 
the literature on numeric methods for two-phase flow [71]-[76] reinforces the advantages of this 
method and the theory on which it is based. 

3.2.1 Hyperbolic Operator 

As was stated in Chapter 2, both phases are compressible, which produces seven distinct wave 
speeds in the system. Equations 2.23-2.24 can then be factorized into a hyperbolic PDE with 
a source term H{U) from the motion of the interface and three separate ODEs for the separate 
components of the non-conservative source term S{U) [31]. 
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(3.36) 



S{U) 


Lg = Lmh ■ Lyji ■ Lpn (3.37) 

When the first-order Godunov seheme is used, factorization of the two-phase model occurs 
according to Equation 3.39. 

jjn+l ^ (3_3g^ 

When the second-order Godunov scheme is used, factorization of the two-phase model occurs 
according to Equation 3.40 [76]. 

jjn+l ^ j^At/2^At^At/2jjn (3 39 ^ 

The ability of an approximate numeric method to capture the shock discontinuity is based on 
the approximate Riemann fluxes at the inter-cell boundary. The method gives a few choices for 
conservative approximate Riemann fluxes. Equation 3.41 is the analogous Godunov update in 
Equation 3.18 for the two-phase system with a first-order Rusanov flux used to update the mass, 
momentum and energy of the gas and liquid phases. 


= C/“ - 


(3.40) 

Fi-\-i/2 = F {Ui^ Ui+i) 

= i (F. + F,+i-S (C/,+,-C/.)) 

(3.41) 



(3.42) 

dx 

2 ■ Ax 

= a” ■ \V-^- { 

) 1 c A// ) 

~ ^g,i ) + ^i-1/2 [C^g,i - ) 

(3.43) 
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The Godunov method for the two-phase system in Equation 3.41 ean be extended to seeond 
order in an analogous was as that deseribed in Seetion 3.1.3. 

3.2.2 Source Term Operators 

Next comes handling of the source terms. They can be factored into three separate terms based 
on timescale or phenomenon and solved sequentially as three separate ODEs after the hyper¬ 
bolic Godunov solution. Eirst the pressure relaxation operator in Equation 3.38 is isolated and 
the solution is given in Equation 3.44 where is the solution vector after the Godunov inte¬ 
gration. 


JJPR Lpr{U^) (3.44) 

This ODE simplifies to Equation 2.41 and 2.42 which are recalled below. 

(agPgEg)^^ - g P g Eg)^ = / PiCia g (3.45) 

JaO 

(aipiEif^ - (aipiEif = - / ' P^dag (3.46) 

JaO 

Internal energy is adjusting in each phase due to the work done by inter facial pressure. The 
goal is to calculate such that the pressures are equilibrated Pg = Pi for all Xi while also 
satisfying the integral constraints in Equations 3.45 and 3.46. A solution procedure used in [31] 
is summarized here. 

L Make an initial guess for based on the difference Pg — Pi 

2. Compute (a^, pg, Eg)^^ and {ai,pi, Ei)^^ from the integrals in Equations 3.46 and 3.47 

3. During relaxation, total mass (a^, pg) is constant so a new iterate of pg can be computed 
and then Eg and Cg as well. The procedure for liquid is identical. 

4. Calculate Pg and Pi from the equations of state 

5. Repeat until equilibrium condition Pg = Pi is reached for all Xi 
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The velocity relaxation, or drag, operator taking the form of Equation 3.47 is next solved nu¬ 
merically. 


jjVR ^ jjPR ^ (3 47) 

The drag force is based on the empirical drag law for spherical objects as given in Equation 
2.43 and shown in discrete form in Equation 3.48. 


Fd{xi,r) = Cdpg{xi,F)Dl{xi,F)ai{xi,F) {ug{xi,F) - ui{xi,F)f (3.48) 

With the drag force on the droplet known, the solution to the ODE is straightforward, as shown 
in Equation 3.49. 


Ug{Xi,F^^) = Ug{Xi,F) -f 
Ui{Xi,F^^) = Ui{Xi,F) - 


At 


ag{Xi,F)pg{Xi,F) 

At 


Fd{xi,F) 


Eg{xi,F+^) = Eg{xi,F) -f 
Ei{xi,F+^) = Ei{xi,F) - 


ai{xi,F)pi{xi,F 

At 


ag{xi,F)pg{xi,F) 

At 


Fd{xi,F) 

Fd{xi,F)Vi{xi,f 


ai{xi,F)pi{xi,F 


Fd{xi,F)Vi{xi,F 


eg{xi,F+^) = eg{xi,r) + -Ug{xi,F^y 
ei{xi,F+^) = ei{xi,F) - ]^ui{xi,F"^^f 


(3.49) 


The third and final source vector left is due to the mass and heat exchange of the phases due to 
vaporization and interface heat transfer. 


Un+l 


= U^^ + Lmh{U^^) 


(3.50) 
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Vaporization dynamics of the droplets is given by the Empirical-Beta Vaporization Law with 
experimental constants given [7]. The present calculation assumes that droplets are mono- 
dispersed and spherical. Their initial size is denoted 77° (xj, 0). When surrounded by hot gas, 
the droplets will start to vaporize and their diameter will evolve according to Equations 3.51 
and 3.52. 


dD'^ 

^(a;„0 = -/7(T,(a;„r)) 


(3.51) 


(3 {Tg{xi,t^)) = 7600 (l + 7.4 ■ 10“^ {Tg{xi,e) - ?,mK)/s 


(3.52) 


Once the new size of the droplets is known, the added amount of gas volume fraction is known. 


4 / £)n+l-^ Dn3^ 

a^+1 =a^+ -ttI p -^ I iV" 

^ ^ 3 2 2 ' P 


(3.53) 


The critical quantity for sinking energy from the gas via vaporization is m which is given in 
discrete form below based on Equation 2.44. 


m{xi,r^^) = pH20v{Pg{xi,r^^)/atm) (ag{xi,r+^) - ag{xi,r)'^ (3.54) 

The density of water vapor changes an order of magnitude in the pressure range of 1-10 atm. It 
is calculated based on a quadratic fit to steam table data [42]. 


PH20v{Pg/atm){Xi,P^^) 


-.0024344- 


Pg{x^,r^^y 


atm 


+.5368 


Pg{xi,t^^y 

atm 


-.077246 kg/m^ 


(3.55) 

With the additional mass of water vapor in the gas phase, the gas density will increase according 
to Equation 3.56 while the density of liquid water does not change. 


p+ = i» 


(a”/)” + m)/a 


n+1 

9 


(3.56) 
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The time integration can now be performed, completing the numerical solution of the source 
vector S{U) at each time step. 

^ 

(agPgUg)^^ = O^g^^ P^g^^^ ^ T TTlVi 

(agPgEg)^^ = p^g^^^ + m {Ln, + E^) + Atg, 

, .MH n+1 n+1 • (3.57) 

[aipi) =ai^ pi^ -m 

{aipiui)^^ = - mVi 

{aipiEi)^^ = a^+^pfEY^ - m (L,, + E,) - AtQ, 

3.3 Solution Procedure for Optimal Control System 

The method for obtaining the discretized flow field, single- and two-phase, was described in 
Sections 1 and 2 of this Chapter respectively. To proceed with numerical conditions based on 
the optimal control system, it is assumed that an entire flow field Uk{xi,t"') is known for all 
spatial indices i, temporal indices n and variable indices k. The vast majority of the additional 
computational complexity required to satisfy all conditions of the optimal system is solving the 
adjoint PDE. Since the method is nearly identical for both the single- and two-phase calculations 
it will be described first. 

A survey of the literature on numerical methods for optimal control [77]-[89] gave insight on 
how to construct an adjoint-based method of solution. However, there were no methods specif¬ 
ically for distributed control with free initial and final data and final time for unsteady shock 
attenuation. 

3.3.1 Discrete Form of Adjoint System 

The continuous form of the adjoint PDE system is shown as the pairing with 5Uj in Equation 
2.69. The system is linear in the adjoint variables. The time derivative of the entire adjoint 
vector for all space at a discrete time E is shown in column-upon-column form in Equation 
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3.58. Let m be the number of spatial grid points and k be the number of adjoint variables. Then 
the adjoint vector L at a discrete moment in time will be of size km by 1 and the matrices in 
Equation 2.69 will be of size km by km. A single component of the adjoint vector, eg. Vi (x,t), 
will be size m by 1 at each time step. 




/ v^(x,,r) \ 


\ 


Vi(x2,r) 




V Vi(Xm,k^) j 

i ^2(a:i,r) \ 

V2{x2,r) 

V ) 


(3.58) 


Vfc(a;2,r) 


\ V Vk{xm,k^) J J 


All of the matrices in Equation 2.69 have a block diagonal structure. The Jacobian matrix for 
all of space at is shown in Equation 3.59. 


57 



Aku{U^) -H- 
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f Aii{xi,D) 
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f Aik{xi,U) 
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\ 

Aii{Xm,U) ) 


[ 

Aik{Xm:U) ) 




^ Aki{xi,U) 

) 


< Akk{Xl,U) 

) 


V 

[ 

Akl{Xm:U) ) 


[ 


/ 


(3.59) 

The adjoint PDE is given in diserete form with an explieit integration in Equations 3.60 and 


3.6E 




(3.60) 


R{U) = A- DU + 


dA 

dx 


dA dU dS 
dU dx dU 


(3.61) 


The matrix D is made up of diserete spatial derivative bloek matriees, eentral differeneing in 
the domain interior, upwind differeneing at the outlet and downwind at the inlet. A single bloek 
is shown in Equation 3.62. 


DmxM — 


1 

2Aa; 


^ -2 
-1 


2 

0 


V 


1 

-1 0 
2 




1 

- 2 / 


(3.62) 


Eaeh time step of the adjoint solution has four parts. Before time integration, the matrix R must 
be assembled. Some of the matriees whieh make up R are known in elosed form and require no 
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discretized derivative. The 3-component system requires assembling A and from the known 
state data. The 7-component system requires, in addition, assembling which has terms with 
mixed closed form and discrete derivatives. The second part of the solution requires assembling 
the matrix and vectors that have discrete derivatives ^ and These two parts can be 
done in parallel. The third part, calculating R, requires sharing memory between processes 
and does not lend itself well to parallelization. With careful direction of memory, there is more 
potential in the assembly of R for speed optimization than will be shown in the results in Section 
4.4. The final part of the adjoint time step is the time integration which boils down to matrix 
addition and matrix-vector multiplication for the explicit scheme. These operations are known 
to be adaptable to parallelization in a straightforward way. 

For adjoint calculations of a scalar PDF with a discontinuity, it has been shown [91] that a 
relaxed system with second order dissipation will recover the non-linear PDE in the limit of 
vanishing viscosity. A small numerical viscosity can stabilize the adjoint solution. These ideas 
have been extended to fluid dynamics systems [92] and are implemented in the current work, 
in a manner which maintains consistency, for both the single- and two-phase numerical adjoint 
solutions. 

3.3.2 Formulation and Solution of Single-Phase Control Problem 

The goal of the calculation is to minimize the finite approximation to the cost functional J while 
also converging to a solution that satisfies all of the necessary conditions of the optimal system. 
The cost functional from Equation 2.79 is given in discrete form in Equation 3.63. 

J = + E^^.fzQ^^max (^Rg(xi,t^) - Q(xi),oJ Ax (3.63) 

Henceforth, i will replace j as the spatial index. The state vector [/ is constrained to the Euler 
System with a distributed control 2 : in the energy balance equation. 

[/ = (p(xi, t'^),pu(xi, t'^),pE{xi, r))^ (3.64) 
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( P \ ( P"^ \ ( 0^ 

^ PM pu^ + P = 0 (3.65) 

\ pE ) \u{pE + P) ) \ z{xi,r) j 

2 

pE{x,,n = pe{x„r) + (3.66) 

P{xi, t”) = pe{xi, P) (7 - 1) (3.67) 

Initial Conditions are stationary air at sea-level throughout the domain. 

p{xi,P) = Ikg/rn? 

u{xi,P) = Om/s 

P{xi,P) = 10^ Pa 

(3.68) 

pu{xi,P) = Qkg/nP-s 

pe{xi,P) = 250000 J 

pE{xi,P) = 250000 J 

The control solution procedure requires a solution of the flow under the influence of the 
control iterate zK Therefore, the initial control iterate will be zero z^ = 0. The first iterate of 
the final time T° is chosen so that the shock wave is allowed enough time to almost traverse 
the domain f2. In this way, the greatest potential for the effect of control action on pressure 
attenuation is possible. At the final time the necessary condition on the adjoint vector was 
given in Equation 2.97 and is shown in discrete form in Equation 3.69 for all Xi G kls- 
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= b ■ max {xi^T^'^ 

, s M (xi, 



= —b ■ max (^(^P (xi, T’’ 

) - Q{xi)) , O) u (xi,T'^) (7 - 1) 

(3.69) 


= b ■ max {J^P {xi^T^^ 

^ — 1 

1 

1 



The adjoint PDE is then solved backward in time as described in Section 3.3.1. This gives all 
of the adjoint variables Vk{xi, t”) over time on the same discrete grid as the flow variables. A 
new control iterate is obtained from the discrete form of Equation 2.99. 

6z^+\xi,r) = -Vi{xi,t^)/a 

z''~^^{xi,t'^) = z^{xi,t^) + 6z^~^^{xi,t^) 

Physical control constraints are then imposed. Since the calculation concerns shock attenuation, 
extracting energy from the gas using a sink is of interest. Consequently, z^ <0 for all space and 
time is enforced. It has been shown that the adjoint variables will travel along the characteristics 
of the flow in the opposite direction. This means that the calculation will suggest putting control 
ahead of the shock wave which we are not interested in since droplet-shock interactions only 
take energy out of the gas behind the shock wave. Therefore, Equation 3.72 gives the restrictions 
on the control. 


Z^^^{Xi,f' 


0 if z'-~^^{xi,t'^) > 0 

0 if Xi> S 


(3.71) 


Easily, the final time is updated according to Equation 3.72 and 3.73. It is important to mention 
that changing the final time does not change the number of time steps N. This means that the 
time step will be slightly different between solution iterates but will still obey the CEE condition. 
This also means that control strengths defined on the time discretization of a previous iterate 
will be shifted when applied to the current flow solution. Overall convergence of the algorithm 
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ensures that these differences are not detrimental and, as the final time converges, so will these 
discrepancies. 
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b{Pg{xi,t^) - Q{xi)) 


(7 - 1 ) ■; 


u, 


Ug (Xi, Pg (Xi, - Pg (xi, 

2 At 

/ PgUg (xi, - PgUg (xi, 

[Xi^r ) -^^^- ^4 


7V-1\ \ 


At 


PgEg (xi, t^) - PgEg [xi, ^ 

At 

Pg (Xi,t^^ - Pg (xi, t^“^) 

At 


]+ 


Ax 


(3.72) 


rj^l + l _ rj^i 


, s{x) 

df , 


(3.73) 


As long as each iterate of the final time is not too far from the optimal value T*, the 
variable will converge in the optimal control solution procedure and the properties in Equation 
3.74 will be true. The definition of % and its discretization are given in Appendix C. 


lim;^„oT'+i = T* 
lim,^o,/(T'+i) = /(T*) = 0 


The discrete form of the functional / from Equation 2.98 is given in Equation 3.72, using the 
definition of the derivatives of pressure from Equation 2.97. The overall algorithm which is 
used to satisfy all of the above necessary optimal conditions is shown in block-diagram form in 
Eigure 3.3 
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Blast wave simulation data at a fixed point in space will be provided as an inlet boundary 
condition. If the state of the flow is supersonic, the condition will be explicitly applied. If the 
flow is subsonic, a non-reflecting boundary condition of the outgoing u — c characteristic wave 
will be enforced [93], [94]. 
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Figure 3.3: Block diagram of solution procedure for single-phase control calculation 
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3.3.3 Formulation and Solution of Two-Phase Control Problem 


The solution procedure for the two-phase control problem is very similar to that of the single¬ 
phase control problem. The major distinction is that an adjoint PDE must be solved backward 
in time to fix an initial condition on a free state variable rather than to calculate a distributed 
control over time. The cost functional J, given in Equation 2.100 and shown in discrete form 
in Equation 3.76, will be minimized while all constraints of the two-phase model are obeyed. 


J = (Pg{xi,t^) - Q(a;i), o)^ Ax (3.75) 

There is no running penalty necessary since the cost associated with control action is penalized 
at the initial time and the cost associated with missing the target pressure profile need only be as¬ 
sessed at the final time. Recall the state vector in primitive form, U = (a^, Ug, Pg, pi,ui, Pi)^ 
and the constraint on the free initial data ag + ai = 1 and 0 < a; < 1, throughout the simulation 
domain fl. The rest of the state variables have fixed initial conditions that are given in Equation 
2.104. The dynamical constraints of the two-phase system are solved in their conservative form 
as in Equation 2.24, but since it is convenient to define the target state in terms of the pres¬ 
sure, and the Jacobian has already been provided [31], the resulting necessary conditions are 
simplified. The adjoint PDE is derived from the Jacobian form in the primitive basis shown in 
Equation 2.51. The non-trivial elements of the matrices are given in Appendices A and B. 

Analogous to Equation 3.69 the necessary condition on the adjoint vector at the final time is 
given for all Xj G 12^. 



b ■ max (^(^P (xi,T''^ — Q{xi)^ ,0^ k = 4 
0 


(3.76) 


With the final time condition, the adjoint PDE can be solved backward in time to arrive at a 
value for Vj{xi, 0). Erom the necessary condition given in Equation 2.107, the volume fraction 
of liquid at each point in space will be iterated, within the constraints, via a Newton method 
shown in Equation 3.77. 
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(3.77) 


To solve for the optimal final time, the Transversality Condition for the two-phase system, given 
in Equation 2.115, is define as a continuous function /(T) with the final time as the independent 
variable. Then f{T*) = 0 at the optimal final time T*. As in the single-phase case, T* can be 
solved iteratively with the Newton method from Equation 3.72 using the discretization of / in 
Equation 3.78. 


f(T') = (P, - Q{xt)) - 


(3.78) 


The definition of ^ and its discretization are also given in Appendix C. The overall algorithm 
which is used to satisfy all of the above necessary optimal conditions is shown in block-diagram 
form in Eigure 3.4. 
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Figure 3.4: Block diagram of solution procedure for two-phase control calculation 
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CHAPTER 4: 
Results and Discussion 


4.1 lOP Simulations in 2D using FASTRAN 





Figure 4.1: Simulated Shuttle lOP: Mach number at 1.2 ms 



3.358E-02I 


Figure 4.2: Simulated Shuttle lOP: Mach number at 4 ms 



3.35SE-02> 


Figure 4.3: Simulated Shuttle lOP: Mach number at 10 ms 

Data on the Shuttle’s grain and chamber pressure [95] was given as input to Cequel, [96] a code 
for steady state rocket property calculations. The output gives the temperature and gas velocity 
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at the nozzle exit plane for a given pressure ratio. Initially, ambient eonditions inside the domain 
are present. The ignition sequenee was simulated in 2D using the ESTFastran Solver [97]. 



Figure 4.4: Simulated Shuttle lOP: pressure at 1.2 ms 



Figure 4.5: Simulated Shuttle lOP: pressure at 4 ms 



2.635E+004 


Figure 4.6: Simulated Shuttle lOP: pressure at 10 ms 

Constant mass-flow boundary conditions equivalent to the steady state exit plane of the rocket 
nozzle on the shuttle were used in the bottom center of the domain on the right face of the step 
as shown in Figure 4.1. Mach number is depicted in three snapshots in Figures 4.1, 4.2 and 4.3 
and pressure in Figures 4.4, 4.5, and 4.6. The last frame is roughly 10 ms after ignition. The 
bottom left edge of the domain represents the rocket body while the bottom right edge is the 
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centerline of the normal to the nozzle exit plane and a symmetry boundary. All other edges are 
non-refleeting boundaries. 

Flow eonditions over time were reeorded at two loeations marked in Figure 4.1. Point 1 is 
near the roeket body 2.5 meters above the nozzle and Point 2 is 1.5 meters along the symmetry 
boundary and the plane normal to the nozzle exit. The eonditions at the recorded locations are 
used as the boundary eonditions in both the single- and two-phase eontrol ealeulations. 

Figure 4.7 shows the flow eonditions over time at Monitor Point 1 near the roeket body and 
Figure 4.8 shows the flow eonditions over time for Monitor Point 2 direetly downstream of the 
nozzle. 
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Figure 4.7: Flow conditions over time for Monitor Point 1 


It is worth mentioning that the two-dimensional flow cannot be truly replieated by an inlet 
boundary eondition to a single-dimensional ealeulation. Negleeting motion in the transverse di- 
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Figure 4.8: Flow conditions over time for Monitor Point 2 

rection while taking (p, u, P) explicitly means that the resulting driver gas has a lower temper¬ 
ature. This is not discouraging since intuition would suggest that hotter gas has more potential 
to vaporize water droplets and hence by the nature of the control action, there is more potential 
for effectiveness. The goal is to be able to handle a range of blast waves that will plausibly be 
seen in an lOP environment. 
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4.2 Single-Phase Control Calculation Results 



Space (m) Time(s) 

Figure 4.9: Monitor Point 2, Pressure profile at final time, no eontrol used. 

These results were published by the author [98]. The solution procedure in Figure 3.3 assumes 
a given target state Q{x) at the final time T*. To illustrate trends in the optimal control solution, 
prescribing a consistent and meaningful target state or sequence of target states is a necessity. 
Figure 4.9 shows the uncontrolled pressure profile over time. An example of a target state, 
shown in red in Figure 4.10, has 85% of the amplitude of the final pressure profile when no 
control is used, OPq, shown in blue. This attenuated pressure profile replaces the uncontrolled 
pressure profile from just behind the shock front to 30 cm from the inlet boundary. The tar¬ 
get is then linearly increased to the magnitude of the uncontrolled pressure profile over 10cm. 
Upstream of the shock, the target state is such that only pressures greater than the target pres¬ 
sure are penalized. This will define in the final time penalty. The nonlinear nature of shock 
waves means that not all target states may be possible for the given initial conditions and bound¬ 
ary conditions. In addition, the purpose of the calculation is to calculate control solutions which 
decrease overpressure and therefore the cost functional should not penalize pressures below the 
target pressure. Furthermore, error near the shock front is not directly penalized but is still 
minimized via the iterative update to the final time. This avoids taking a variation of the state 
variables near discontinuities, which will violate the assumption of a small variation. In each 
of the given examples the weighting constants a and b are 10“® and 10^ respectively. Since b is 
much larger than a, minimizing J is dominated by minimizing the final time penalty. The value 
for b is such that the magnitude of the third adjoint variable V 3 {x, t), and thus the control action 
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z{x, t), are sufficient to cause noticeable attenuation to the blast wave. 



Space (m) Time(s) 


Figure 4.10: Monitor Point 2, Optimal pressure profile at final time, target state Q{x) ~ 
85% OPo 



1.5 

Space(m) 


2 


Figure 4.11: Monitor Point 2, Optimal distributed control,Q(a:) ~ 85% OPq 

The desired state cannot be directly set as described but must be gradually scaled down to the 
desired overpressure to assure convergence of the algorithm. It was found that when starting 
with an overpressure of 8:1 from Monitor Point 2 data, setting a target state with 85% of the 
uncontrolled overpressure, as in Figures 4.10, was the most attenuation to the shock that could 
be desired and still maintain convergence of the the algorithm. The reason for this can be 
understood in the nature of the optimal system. The optimal conditions are all coupled and 
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Figure 4.12: Monitor Point 2, Optimal pressure at final time, Q(x) ~ 70% OFq 
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Figure 4.13: Monitor Point 2, Optimal distributed control, Q(x) ~ 70% OFq 


inter-dependent such that a calculated blast wave under no control action is in no way close to 
the optimal pressure profile of a blast wave under considerable damping control action so as to 
render a blast wave that matches the target state with a significantly diminished overpressure. 
When the final time solution is not nearly optimal, then first-order corrections to these terms 
are not enough to exploit in an iterative sense toward finding an optimal control satisfying all 
necessary conditions. 

Figure 4.11. is the corresponding distributed optimal control solution. Figures 4.12 and 4.13 
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are the optimal pressure and control solutions respectively for a target state that has 70% of the 
uncontrolled overpressure. Figures 4.14 and 4.15 are the optimal pressure and control solutions 
respectively for a target state that has 50% of the uncontrolled overpressure. 
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Figure 4.14: Monitor Point 2, Optimal pressure at final time, Q{x) ~ 50% OPq 
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Figure 4.15: Monitor Point 2, Optimal distributed control, Q{x) ~ 50% OPq 


The top plot in Figure 4.16 is the logarithm of the cost functional J over solution iterations. 
The jumps in J every 25th iteration indicate that the target state Q{x) has been redefined with 
a lower overpressure. The bottom plot in Figure 4.16 shows the corresponding final time T* 
iterates with T° = 2 ms. 
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A negative control in the energy equation can be thought of as an internal energy sink. Tem¬ 
perature is directly related to internal energy, so the control action cools the gas. A cooler gas 
has a lower speed of sound. Consequently, it will take longer for the shock front in the wave 
interacting with control to reach the same point in space that the shock without control reached. 
Additionally, from a hyperbolic wave theory point-of-view, a non-linear wave travels slower if 
it has less amplitude. The control action’s purpose is to attenuate amplitude and this necessar¬ 
ily slows the wave down. The Transversality Condition is very sensitive to error at the shock 
front since this is where the time rate of change of the pressure is greatest. It would be a diffi¬ 
cult matter to get significant attenuation in a converged solution with a fixed final time for this 
reason. 




Figure 4.16: Monitor Point 2, Cost functional (top) and final time (bottom) vs. iteration, Q(x) ~ 
50% OFo 

In each example, a physical restriction on the control is imposed such that energy can only be 
taken out of the gas behind the shock wave. This more accurately portrays how discrete droplet 
sprays will sink energy from the gas via vaporization. The plots shown in Figures 4.17 and 4.18 
are distributions of the energy equivalent vaporized water mass to the optimal control solutions. 
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By dimensional analysis, inspection of the energy balance equation indicates that the units of 
z*{x, t) are Watts. Integrating the control solution from [0, T*] gives an energy distribution in 
space. This energy can be directly equated to an amount of water mass that must be vaporized 
using the latent heat of vaporization of water Lh^ = —2.26 ■ 10® J/kg at 100° C. Equation 4.1 
relates the optimal control solution, distributed in space and time, to an equivalent distribution 
of water mass vaporized over the same time interval. 

1 rT* 

= -— / Z*{x,t)dt (4.1) 

J^hv 



Figure 4.17: Monitor Point 2, mH 20 vape{x), Energy equivalent distribution of water mass vapor¬ 
ized 
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Figure 4.18: Monitor Point 1, mH 20 vape{x), Energy equivalent distribution of water mass vapor¬ 
ized 
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Figures 4.17 and 4.18 show how the equivalent vaporized water mass distribution increases with 
increasing shock front attenuation. It can clearly be seen that more water is required to attenuate 
the stronger shock front using Monitor Point 2 data from directly in front of the rocket nozzle 
than from the weaker shock, Monitor Point 1 data, taken near the rocket body. Water takes time 
to vaporize and therefore cannot sink energy from the gas instantaneously and arbitrarily close 
to the shock front. In practice, water will take energy out of the gas at least a meter behind 
the shock front. That less energetic gas will then drive the shock with less energy and the 
pressure profile will flatten out to a diminished overpressure. With more detailed physically- 
based restrictions on the control action, this single-phase shock wave control formulation can 
be made to better replicate a two-phase shock-droplet interaction control calculation without 
the considerable added complexity of a two-phase compressible flow model. Regardless, trends 
in the data seem heuristically correct and therefore empirical scaling laws in magnitude and 
location may be sufficient. 

4.3 Two-Phase Calculation Results 

The results for the two-phase control problem were obtained with the algorithm shown in Figure 
3.4. Figure 4.19 shows the optimal initial condition for the free water volume fraction variable. 
The curves give the optimal water volume fraction distribution for a given target state. The 
legend tells what percentage of the absolute overpressure of the uncontrolled blast wave was 
used to define the target state. In this case there is no artificial upper constraint on the maximum 
value of the volume fraction, aside from the physical limit of 1. 

The limit of controllability over two meters is close to the target state 76% OPq. For each 
optimal initial water distribution, the resulting pressure profiles at the optimal final time T* are 
shown in Figure 4.20. 
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Figure 4.19: Optimal water volume fraetion distributions af (x, 0), uneonstrained at initial mo¬ 
ment in time 
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Figure 4.20: Optimal final time pressure profiles resulting from initial data in Figure 4.19 
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Figure 4.21: Cost functional J in red and initial penalty I over iterations of solution procedure 
in Figure 3.4 
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Figure 4.22: Final time, in seconds, vs iteration number for the unconstrained case. The regions 
of greatest slope are after iterations where the target state has been changed. The slope tending 
toward zero means that the solution is converging. 
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The convergence of the algorithm is shown in Figure 4.21. Every 25 iterations, the target state’s 
amplitude is ramped down so that optimal solutions that yield ever increasing attenuation levels 
can be determined in a single execution of code. The history of the final time variable over the 
iterative solution procedure is shown in Figure 4.22. 



Figure 4.23: Optimal initial water volume fraction distributions, ai{x, 0) < .7 


Another case was run with an upper bound on the volume fraction of water at 70%. The calcu¬ 
lated optimal water distributions are shown in Figure 4.22 with the resulting pressure profiles 
shown in Figure 4.23. 
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Figure 4.24: The optimal pressure profiles at the final time resulting from the eontrol initial data 
in Figure 4.23 
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Figure 4.25: Optimal initial water volume fraetion distributions, ai{x, 0) < .5 

A final ease was run, this time with an upper bound on the volume fraetion of water at 50%.The 
calculated optimal water distributions are shown in Figure 4.25 with the resulting pressure pro¬ 
files shown in Figure 4.26. 
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Figure 4.26: Optimal final time pressure profiles resulting from initial data in Figure 4.25 
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A few more isolated results were obtained in the interest of categorizing the significance of 
model parameters. Figure 4.27 shows the effect of droplet size on overpressure attenuation. 
Each pressure profile is plotted after two milliseconds of shock propagation from the left bound¬ 
ary toward the right. For a constant amount of water mass, more surface area of the droplets 
are exposed to the flow when the droplets are smaller. This is verified in Equation 4.2 for 
mono-dispersed droplets with diameter Dp. 


S'„ = 47r ( ^ ) ■ N, 


Np ft;- 


2 

Ax^ 




i^p 


(4.2) 


This result agrees with empirical results that more exposed droplet surface area, the greater the 
overpressure attenuation [12]. Since surface area is maximized with infinitely small droplets, 
initial droplet size isn’t an optimizeable parameter. Droplets with a diameter Dp = 25/im were 
used in all of the preceding results since that size is small enough to have a significant effect 
over two meters, is a reasonable size based on injector atomizer specifications and prevents the 
droplets from being completely vaporized within the simulation time T*. Avoiding complete 
vaporization is desirable since a non-zero volume fraction of either phase is a prerequisite to the 
model being used [31]. 
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Figure 4.27: Effect of variable droplet size 
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Water vapor density was shown to have a drastie effeet on results for optimal water distributions. 
Aeeording to the steam table data [42], water vapor density inereases by an order of magnitude 
from ambient pressure to the 8 atm level behind the shoek front. Therefore, if water vaporizes 
under high pressures, mueh more mass ehanges phase and eonsequently, the dissipative effeet 
of vaporization is mueh more signifieant. The results in Figure 4.28 show that more than twiee 
as much water is needed to get 5% attenuation in pressure depending on what the gas pressure 
surrounding the droplets when they vaporize. 



Figure 4.28: Effect of variable water vapor density due to gas pressure 
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Figure 4.29: Effect of maximum of Gaussian water volume fraction distribution on lOP strength 
after 1 ms. 

Figure 4.29 shows that for an initial water volume fraction distribution with a Gaussian shape, as 
the water volume fraction at the inlet approaches 1, the effect of plume confinement [29] takes 
over and can actually increase the jump in pressure at the shock front. The Gaussian distribution 
is such that the maximum exists at the inlet boundary and about 1% of the maximum half-way 
into the domain. 
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4.4 Serial and Parallel Code Run-Time Optimization 


The last section of the results shows the run-time optimization performed within Matlab. For 
both the single-phase (3 adjoint) and the two-phase (7 adjoint) systems, 5 cases of increasing 
size were carried out with serial and parallel syntaxes. Jacket [99] was used to run the 3 adjoint 
calculation in parallel on Nvidia GPUs. Results for the time duration needed to solve the system 
in [T*, 0] are shown in Figure 4.30 and those of the 7 adjoint system are shown in Figure 4.31. 


Case 1: 
Case 2: 
Case 3: 
Case 4: 
Case 5: 


m = 
m = 
m = 
m = 
m = 


100, timesteps = 1000 
200, timesteps = 2000 
300, timesteps = 3000 
400, timesteps = 4000 
500, timesteps = 5000 



10^ discrete points in space-time 
4-10® discrete points in space-time 
9-10® discrete points in space-time 
1.6-10® discrete points in space-time 
2.5-10® discrete points in space-time 



Figure 4.30: Time duration for solving the 3-component adjoint system of the Euler equations 
over the time interval [T*, 0] 
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Figure 4.31: Time duration for solving 7-oomponent adjoint system of two-phase model over 
the time interval [T *, 0] 
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CHAPTER 5: 
Conclusions 


A new iterative solution procedure was developed, which can calculate optimal distributed con¬ 
trol solutions for systems of quasi-linear hyperbolic partial differential equations with free initial 
and final data and final time. This procedure has been successfully applied to single- and two- 
phase compressible gas dynamics in one dimension with the goal of diminishing overpressure 
at the shock front of a blast wave generated by an ignition overpressure. Examples of opti¬ 
mal attenuation to blast waves typically encountered in the launch environment of the Shuttle’s 
SRBs during an ignition are given. The single-phase control solutions can be seen as a spatial 
distribution of energy equivalent vaporized water mass required to achieve a given level of over¬ 
pressure reduction. Results for these mass distributions are shown for inlet boundary conditions 
like those given in the monitor point data. 

The same methodology for solving the control system is implemented on a two-phase model. 
The control action takes the form of the free initial data describing a distribution of water 
droplets. Optimal water volume fractions are calculated for increasing levels of attenuation. 
Cases where the maximum volume fraction of water is restricted to 50% and 70% are also 
presented. 

The results give several key insights relevant to implementing water injection systems for blast 
wave attenuation. Smaller droplets will vaporize quicker since the total surface area exposed 
to the flow is larger [12], resulting in cooler gas driving a shock with an attenuated jump at 
the shock front. This would suggest that large regions of space completely filled with water 
are sub-optimal for blast attenuation since connected streams or ligaments of water expose less 
surface area to the gas. 

The degree of attenuation depends largely on the rate of mass changing phase from liquid water 
to gaseous water vapor due to forced vaporization m, which is equal to the product of the water 
vapor density and the rate of change of the volume occupied by the water droplets. In the high 
pressure and high temperature region behind the shock front, water vapor density is high and the 
rate of vaporization is high as well which both contribute to a large value for m. As the control 
takes effect, the pressure and temperature of the gas both decrease meaning that the effects are 
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being felt by the shock front at a slower rate. These results show that, in general, stronger 
shocks can be attenuated more rapidly than weaker shocks via water droplet vaporization. 

Similar to work that has identified plume confinement as being disadvantageous to blast attenu¬ 
ation [29], the results also show that when a shock propagates through a water volume fraction 
very close to 1 over a sufficient interval in space, so much liquid water is converted to gaseous 
water vapor that the increase in gas density becomes more significant to the pressure than the 
cooling of the gas due to vaporization. This shows that for a given interval of space there is 
a limit to the controllability. In a two meter domain, 76% OPq was the target with the low¬ 
est overpressure for which the solution still converged. The case with the volume fraction of 
water restricted to a maximum of ai < 70% converged with the lowest target overpressure 
of 80% OPq. With the water restricted to below ai < 50% by volume, the most attenuation 
achievable with a two meter domain is 82.5% OPq. 

The code developed by the author is still in the prototype stage, however, the run-time opti¬ 
mization results show that certain parts of the code have been optimized within Matlab script¬ 
ing. Implementation of the two-phase calculation on parallel GPUs using Jacket or through 
other GPU platforms would make the code a feasible tool to use in lOP and other two-phase 
unsteady shock wave attenuation calculations. The framework of the control formulation and 
method of solution will remain unchanged while allowing for added degrees of freedom like 
poly-dispersed droplets, droplet breakup, surface tension and multiple dimensions. 
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APPENDIX A: 

Calculation of the elements of for the 

two-phase adjoint system 


Row 1: 


Row 2: 


Row 3: 


9 A 
^ 

OUk 

dx 


^ dVjdpi ^ dVjdui 
dag dx dpg dx dug dx dpi dx dui dx 


(A.l) 


d _ dtik 
duk dx 


Pg /dv^ _ Vj - Ug /dp^ 

Ug \dx dx J ag \ dx 


p^dog 
ag dx 


(A.2) 


d duk _ dug 

X. - ^22 ■ — 

UU}^ ux ux 


(A.3) 


d duk _ dpg 

7^-^23 * 

UUk UX UX 


(A.4) 


d _ duk 
duk dx 


-1 ^ + (p - Pi) f—— + —— 

agPg dx dx ^ * V'^C/ Pg 


(A.5) 


d duk dug 

^—A33 ■ —— — —— 
uU}^ UX ux 


(A.6) 


d _ ^ ^ -l<9pg 
duk dx pI dx 


(A.7) 
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Row 4: 


•9 ^ du]^ Pg^gi f dVi d)Ug Vi Ug ( 9 q ;^ \ ^ Vi Ug /, dPi ^ d)Pg 

duk dx Ug y dx dx Ug dx ) ag \ dx dx 


d duk dPg 

7^^43 ■ = 777^ 

UUj^ <JX uX 

d . duk _ 

duk dx dx 

Row 5: 

d duk _ Pi (dVj dui \ ^ Vj-ui (dpi ^ pi dag 

duk dx 1 — ag \ dx dx J 1 — ag \dx 1 — ag dx 

d duk _ dui 

xi .—^55 ■ — 7^ 

UU]^ ox ox 

d duk _ dpi 

-7s -^56 ■ 7^ — 7^ 

uUk dx ox 


Row 6: 


d . _ d^ 

duk dx 


-1 

(1 - aj Pi 


(m 

dx 


m 

dx 


+ Xi - p^) 


-1 dag ^ 1 (9p/\ 

1 — ag dx Pi dx J 


d duk _ dui 

Xs - ^66 ■ 7 ^ — 7 ^ 

uUk OX OX 

d _ ^ ^ -I dpi 
duk dx pf dx 


(A.8) 

(A.9) 

(A.IO) 

(A.ll) 

(A.12) 

(A.13) 

(A. 14) 

(A.15) 

(A.16) 
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Row 7: 


duk _ Picl (DV, dui V,-Ui dag\ V,-Ui f DP, dPA 

duk dx 1 — ctg y Dx dx 1 — ag dx J 1 — y Dx dx J 

(A.17) 


d duk _ dPi 

—A76 ■ — 7 z“^ 

UUj^ UX ux 


(A.18) 


d duk _ ^ 

duk dx dx 


(A. 19) 


109 



THIS PAGE INTENTIONALLY LEET BLANK 


no 



APPENDIX B: 

Calculation of the elements of for the 
two-phase adjoint system 

dS 

Solution of the adjoint system of PDEs requires deriving the matrix where S (U) is the 
source vector for interactions between the two phases. 


/i {Pg - Pi) + Aag 


mVi + Fd 

S{U) = m {Lhy -f Ei) -\- Qi -\- FdYi — ^-Pi {Pg — Pi) 


(B.l) 


-mVi - Fd 

^ -m {Lhv + Ei) - Qi- FdVi + fiPi {Pg - Pi) I 

Before calculating derivatives with respect to state variables, it must be made clear which terms 
are explicit functions of which variables. Aagyap {U ), m {U ), Fd {U ), Qi {U ), Vi {U ), Pj {U ), 
and Ei{U) yet only Ei is an explicit function of all seven state variables. Therefore, it is worth 

determining only the derivatives that exist ahead of time to more precisely define the elements 

,dS 

The easiest term to differentiate is the interface pressure. 


Pi — OigPg + (1 ~ (^g) Pi 


(B.2) 


Therefore, Pi = Pi (a^, Pg, Pi) and the derivatives are: 


dP 


(B.3) 
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(B.4) 


dPi 

dP, 


= a 


9 


dP^ 

dPi 


1 ag 


(B.5) 


The interface velocity is a function of a subset of state variables as well: Vi = Vi (a^, pg,Ug, pi,ui) 
The derivatives are: 


dag 


PgUg - PlUl 

^gPg T (1 Pi 


^gPg^g *^ 9 ) Pi^i / \ 

{o^gPg + {I - o^g) Plf ' 


(B.6) 


dVj _ OigUg _ O^gPgUg + (1 - Pig) PlUl ^ 

dpg agPg + {I - ag) Pi (+ (1 _ «J ^ ^ 


(B.7) 


dVj _ (1 - ag) Ui CtgPgUg + (1 - Olg) PlUl 

dpi agPg + {1 - ag) Pi [agPg + {1 - ag) Pl)^ ^■ 


(B.8) 


dVi _ agPg 

dug agPg T (l ^9) Pi 


(B.9) 


dVi _ (1 ^g) Pg 

dui agPg + (1 - ag) Pi 


(B.IO) 


The drag force Fd iU) = Fd (a^, pg, Ug, m) 


d_Fd 

dag 


3 

4 


^dPgPp i.^g 


Ul)^ 


(B.ll) 


m 

dpg 



(^g) K - '^if 


(B.12) 


dPd 

dUg 


\c,f>,Dl (1 


ag) {Ug Ul) 


(B.13) 
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3 

2^d,P9^p ~ ^ 9 ) (^9 ~ 


(B.14) 


dFd 

dui 


The derivatives of m and Aag^ap are elosely related. Reeall, 


rh iu) = Attgyap iu) ■ PH 20 AP 9 ) (B.15) 

where is the density of water vapor and is only a funetion of gas pressure. It’s derivative 
is easily ealeulated from the quadratie fit to the empirieal steam table data given in Equation 
2.38. 


^ f.5368 - .0048688 ■ kg/m^ (B.16) 

dPg / atm \ atm) 

It suffiees to ealeulate derivatives of Aagyap with respeet to the state variables. The vaporization 
model is based on the Empirieal-Beta Vaporization Eaw whieh states: 


(d”+‘)" = {pi)‘ - A0 (r;) 


(B.17) 


and 


^ ((C”)' - ((C”)' - Atl3 (p„ (B.18) 

Tg = where R is the speeifie gas eonstant for air, so Aagyap {U) = Aagyap {Pg, Pg) ■ The 
derivatives are: 


dAa 


gvap 


dpg 


NpTT /3\ 
6Aa;3 V2/ 




1/2 


(At) 


dp 

dpg 


(B.19) 


and similarly, 
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dAagyap _ NpTT f3\ ({ f \J.\ 


dP„ 6Ax^ V2 


{d;) -At/3) (At) 


dPr, 


(B.20) 


The vaporization term is given by Equation B.21: 


13 {pg, Pg) = 5o + Bi - Tmi^ j pni^/s 


(B.21) 


where the constants are given as Bq = 7600, Bi = 7.4 ■ 10 B 2 = 2.7548 and T^m = 300 
degrees Kelvin. Therefore, the derivatives are: 


T-> l-> l-> f rn rn \B 2 — 1 ('^9 


dpc 


— —BqBiB 2 {Tg — T„ 


.Pg. 


(B.22) 


TJ T2 T2 /rn rn \ ^2 — 1 (^9 


dPy 


— BqBiB 2 {Tg — T„ 


Pr, 


(B.23) 


The total energy at the inerface is simply defined in the conservative basis: 


Bi — OigEg + (1 — ag) El (B.24) 

But not so easily in the primitive basis. Using the equations of state: 

= “» - “»> - ■®“‘) 

In this basis, the interface velocity is a function of all primitive variables and therefore there are 
seven non-zero derivatives to calculate: 


dag 


Eg - El 


(B.26) 
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(B.27) 


dEj _ -agPg 

dpa pI (79-1) 


dEj _ - (1 - ttg) {Pi - Jim) 

dpi pf (7z - 1) 


(B.28) 


dUg 


— CHgUg 


(B.29) 


dui 


(1 - ttg) Ul 


(B.30) 


dEj _ Ug 

dPg Pa (79 - 1 ) 


(B.31) 


dEj _ I - ag 

dPi pi ( 7 / - 1 ) 


(B.32) 


Lastly, the derivatives of the rate of eonveetive heat transfer between the two phases at the 
interphase Qi (U) are ealeulated. 


Q, {U) = h-Sg- {Tg - Tl) 

= ^Sg-{Tg-Tl) 

= (2 + .459PrV3i?e-^5)^.49r(t)'-(5, 


1 ( Pl+llT^l 

Cvi V Pi (71-1) 



= 2 + 


450 ^ DpPgjUg—Ul) 



1 / Pl+'ll'^l 

Cvi VPi(7i-1) 



Therefore Qi = Qi {pg, Ug,Pg,pi,Pi) and the derivatives are: 


(B.33) 


dQi 

dpa 


XnDp 



Tl) - Nu^\ 
Pa) 


(B.34) 
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(B.35) 


dQi 

dUg 


XnDp 


dNu 

dug 


(Tg 


Ti) 


dQi 

dPg 


XnDpNu 


P 9 


(B.36) 


dQi 

dpi 


T 

XttDpNuP- 

Pi 


(B.37) 


dQi 

dui 


XnDp 


dNu 

dui 


(Tg 


Ti) 


(B.38) 


dQi 

m 


XjrDpNu— 

Pi 


Ti _ 

- llT^l 


And it can easily be seen that the derivatives of the Nusslets’ Number are: 


(B.39) 


^ = .458 ■ .55 ■ ^i) 

dPg P 


(B.40) 


^ = .458 ■ .55 ■ (B.41) 

OUg P 


^ = -.458 ■ .55 ■ Pr-P^Re-^^-^^ (B.42) 

oui p 

Now, with all of the derivatives of eaeh eaeh individual funetion of the state veetor that eompose 
the souree veetor, the on-zero elements of the matrix ean be ealeulated. 

Row 1: 


dSi 

dpg 


NpTT /-3\ 
6Aa;3 V 2 / 




1/2 


d(3^ 

dpg 


(B.43) 
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(B.44) 


Row 2: 


Row 3: 



dag dag dag 

dSs ^ ^ I • ^ I ^ 
(9pg dpg * ^(9pg (9pg 

9^3 _ . < 91 ", 

dpi "^dpi 
^ = rh—+ — 

dUg dUg dUg 

^=rn— + — 
dui dui dui 


(B.45) 

(B.46) 

(B.47) 

(B.48) 

(B.49) 

(B.50) 

(B.51) 

(B.52) 


dSs dm 
Wg = Wg^^ 


(B.53) 
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Row 4: 


dS4 dEi dF., ^ dVi 

Oag OUg OUg Ottg 


dSi _ dm , ^dPi , dQi , dPd^, , ^ dVi 

a— “ a— \^hv + Ei) + m— -h —-h ——Vi + Fd -^— 

^Pg ^Pg ^Pg ^Pg ^Pg ^Pg 


dS, _ .dE. dQ, , dFd^^ , ^ dV 

dUg dUg dUg dUg dUg 


dSi _ dm ^ ^ ^dEi ^ dQi ^ D^ , D^ 

dPg dPg^^^^'^^^^^'^dPg^ dPg ) 


dS, ^ .dE, dQ, dV 

dpi dpi dpi ^dpi 


dS, . dE, dFd ^ ^dV. dQ, 
= m— + Vi— + Ed— + 


dui dui dui dui dui 


dSi . dEi dQi 


(B.54) 

(B.55) 

(B.56) 

(B.57) 

(B.58) 

(B.59) 

(B.60) 


Rows 5, 6 and 7 are the opposite sign of rows 2, 3 and 4 respectively. 
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APPENDIX C: 

Time Derivatives and Discretization of 
Transversality Condition for Single- and 
Two-Phase Control Formulation 

For the single-phase problem, the eontinous form of the funetional / that must equal zero at 
the optimal final time is given in Equation 2.98. For eonvenianee, the following definition is 
eonstrueted, suppressing the subseript g for the gas phase. 

dP dp 
dp dt~'^ 
dP dpu 
dpu dt 
dP dpE 
dpE dt 
dP 

The diserete form is shown in Equation C.2. 

2 At 

u{xi,r) - — -+ 

pE{xi, t^) - pE{xi, t^~^) 

At 






dt 



^dP dP 
dt dt 


+ b{P-Q) 


d^P\ 

dP j 


dx 


(C.3) 


The eontinuous derivative of the term ^ is given in Equation C.4 and the diserete form in 
Equation C.5. 
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dudp u^d‘^p dudpu d‘^pu d‘^pE\ d^P 

dt dt^ 2 dP ^ dt dt ^ dP dP j dP 


= (7 -1) 


u{xi,t^) 


u{Xi,t^) - u{Xi,t^ p{Xi,t^) - p{Xi,t^ 


” ^ At At 

P{xi, t^) p{xi, t^) - 2p{xi, t^~^) + p{xi, t^“^) 

2 AP 

u{xi, t^) - u{xi, pu{xi, t^) - pu{xi, 

At At 


M(a;i,t^) 


^ pu{xi,t^) -2pu{xi,t^ P + pu{xi,t^ 2 ) 


pE{xi,t^) - 2pE{xi,t^ P + pE{xi,t^ 

At2 ^ 

1 P(a;,, t^) - 2P(a;,, t^"^) + Pjxj, t^"^) 

7-1 At 2 


(C.5) 


df 


z{Xi,t^) - z{Xi,t^ P{Xi,t^) - P{Xi,t^ 

At ^ At 

b{P{x„t^)-Q{x„t^))-f^DP 


DP+ 


(C.6) 


For the two-phase problem, the derivative of the eontinuous form of the Transversality Condi¬ 
tion is ealeulated from Equation 2.115 and shown here in Equation C.7 and it’s diserete sum¬ 
mation form in Equation C.8. 


df f fdpy 


(C.7) 
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^ r^yM 

dt ^ 




(^P{xi,t^) - Q(xi)J 


At 


+ 


F(xi, t^) - 2P{xi, + P{xi, t^-2) 

AP 


Ax 


(C.8) 
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